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ABSTRACT 


The  major  objective  of  this  project  is  to  design  simple  low-dispersion  Finite-Difference 
Time-Domain  (FDTD)  methods  for  electromagnetics.  Literature  review  indicated  that  the 
Nonstandard  Finite  Difference  (NSFD)  method  exhibits  great  potentials  in  dispersion  re¬ 
duction.  Different  from  the  Standard  Finite  Difference  (SFD)  methods,  the  NSFD  methods 
are  derived  directly  based  upon  dispersion  analysis.  In  this  report,  the  basic  concepts  of 
the  NSFD  methods  are  generalized  to  various  extended  finite-difference  stencils.  Further¬ 
more,  several  improved  NSFD  methods  are  presented  based  on  the  standard  fourth-order 
stencil  to  mitigate  the  dispersion  in  multi-dimensional  domains.  The  least  square  method 
is  used  to  optimize  the  coefficients  of  the  proposed  stencils.  Numerical  simulations  show 
that  these  schemes  significantly  reduce  the  dispersion  error  of  their  standard  counterparts. 
Many  technical  issues  in  practical  implementations,  such  as  absorbing  boundary  conditions, 
stability  conditions,  and  Gauss’s  Laws,  are  discussed  and  justified.  Moreover,  two  special 
conditions  are  proposed  for  the  extended  stencils  in  the  vicinity  of  the  dielectric  material 
discontinuities.  It  was  demonstrated  that  the  accuracy  of  the  fourth-order  stencil  is  fully 
restored  by  applying  these  conditions. 

A  very  important  electromagnetic  interference  (EMI)  problem  is  also  presented  in 
this  report.  The  problem  concerns  the  effects  of  passengers  on  personal  electronic  device 
(FED)  mutual  coupling  in  a  simplified  Boeing  757  fuselage.  The  results  predicted  by  EDTD 
are  compared  to  the  measurements  performed  in  the  Electro-Magnetic  Anechoic  Chamber 
(EMAC)  of  Arizona  State  University.  It  was  found  that  the  presence  of  the  passengers 
significantly  dampens  the  resonances  in  the  fuselage,  which  lower  the  potential  EMI  threat 
to  the  on-board  communication/navigation  systems. 
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CHAPTER  1 


INTRODUCTION 

Electromagnetic  theory,  one  of  the  most  important  branches  of  physical  science,  stud¬ 
ies  the  electric-magnetic  fields  generated  by  charges  at  rest  or  in  motion  (currents)  [2]. 
The  relations  and  variations  between  the  electric- magnetic  fields,  charges  and  currents  are 
governed  by  basic  laws  developed  by  a  group  of  distinguished  scientists:  Earaday,  Ampere, 
Gauss,  Lenz,  Coulomb,  Gauss,  Volta,  and  others.  In  1873,  Maxwell  was  able  to  mathe¬ 
matically  combine  these  laws  into  a  consistent  set  of  vector  equations,  which  are  referred 
to  as  Maxwell’s  equations.  The  modern  form  of  Maxwell’s  equations  were  formulated  by 
Heaviside  in  the  late  19th  century  [3].  Then  Hertz  completely  validated  the  theory  of 
Maxwell’s  by  a  series  of  experiments.  Motivated  by  the  practical  applications,  most  impor¬ 
tantly  radar  and  wireless  communications,  electromagnetic  theory  has  been  developing  for 
over  100  years.  Today,  electromagnetics  has  become  an  indispensable  discipline  of  everyday 
life  which  covers  a  wide  range  of  applications,  such  as  antennas,  scattering,  microwave  cir¬ 
cuits  and  devices,  radio-frequency  and  optical  communications,  broadcasting,  geosciences 
and  remote  sensing,  radar,  radio  astronomy,  quantum  electronics,  solid  state  circuits  and 
devices,  eletromechanical  energy  conversion  and  computers. 

Almost  all  electromagnetic  phenomena  can  be  well  described  by  Maxwell’s  equations 
either  in  differential  or  integral  forms,  associated  with  appropriate  boundary  conditions. 
However,  only  a  few  canonical  problems  can  be  solved  exactly.  Closed-form  solutions  for 
most  practical  and  complex  electromagnetic  problems  are  usually  intractable.  Eor  example, 
some  problems  involve  complex  geometries  that  do  not  conform  to  any  known  orthogonal 
coordinate  systems  for  which  the  scalar  Helmholtz  equation  can  be  solved  by  the  method  of 
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the  separation  of  variables.  Also,  many  other  problems  may  include  nonlinear,  anisotropic 
or  dispersive  materials.  Consequently,  many  numerical  methods  have  been  developed  to 
solve  such  problems. 

With  today’s  computer  resources,  numerical  methods  are  drastically  more  efficient,  if 
not  the  only  way  to  solve  complex  practical  electromagnetic  problems.  The  most  popular  nu¬ 
merical  methods  include,  high  frequency  asymptotic  methods.  Method  of  Moments  (MoM), 
Finite  Element  Method  (FEM)  and  Finite-Difference  Time-Domain  (FDTD)  method,  etc. 
Brief  introductions  to  these  methods  follows. 

High-frequency  asymptotic  methods  use  optical  asymptotic  techniques,  such  as  Geo¬ 
metrical  Optics  (GO),  Geometrical  Theory  of  Diffraction  (GTD),  Physical  Optics  (PO), 
and  Physical  Theory  of  Diffraction  (PTD).  These  methods  calculate  the  incident,  reflected 
and  diffracted  fields  separately  by  introducing  reflection  and  diffraction  coefficients.  The 
ultimate  solution  is  a  superposition  of  the  contributions  from  each  part.  These  methods  pro¬ 
vide  physical  insight  into  the  radiation  and  scattering  mechanisms  and  yield  accurate  results 
that  compare  well  with  the  experiments.  Pure  high  frequency  methods  are  very  useful  in 
electrically  large  problems  with  simple  materials,  such  as  perfect  electric  conductor(PEG). 
However,  these  methods  are  usually  more  difficult  to  apply  to  complex  geometries,  espe¬ 
cially  those  that  involve  complex  materials.  Thus,  high-frequency  asymptotic  methods  are 
usually  hybridized  with  other  numerical  methods,  such  as  MoM.  An  introduction  to  the 
methods  can  be  found  in  [2]. 

Gompared  to  the  high-frequency  asymptotic  methods,  MoM  can  handle  objects  of  ar¬ 
bitrary  geometries  especially  those  of  planar  shapes.  It  casts  the  solution  for  the  induced 


current  density  in  the  form  of  an  integral  equation,  and  expands  the  unknown  quantity  into 
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a  set  of  basis  functions.  MoM  has  become  popular  since  the  pioneering  work  of  Richmond 
and  Harrington  in  1960s.  Since  then,  it  has  successfully  solved  a  wide  variety  of  prac¬ 
tical  electromagnetic  problems,  such  as  radiation  from  thin- wire  antennas  and  scattering 
problems.  However,  MoM  usually  cannot  easily  deal  with  inhomogeneous  and  non-linear 
materials  because  of  the  difficulty  to  find  the  analytical  Green’s  function.  Details  of  the 
method  can  be  found  in  [4]. 

FEM  originates  from  material  and  structural  analysis  in  mechanics  and  civil  engineer¬ 
ing.  It  was  not  used  to  solve  practical  electromagnetic  problems  until  the  late  1960s.  The 
successfulness  of  FEM  is  due  to  its  powerful  capability  to  model  general  anisotropic,  dis¬ 
persive  and  nonlinear  media.  By  volumetrically  discretizing  all  space  and  transforming 
Maxwell’s  equations  into  a  weak  form  of  either  the  electric  or  magnetic  field  wave  equation, 
FEM  is  able  to  solve  more  complex  problems.  The  mathematical  advantages  of  using  the 
weak  form  of  the  wave  equation  are  that  its  solutions  do  not  have  as  restricted  properties 
as  the  solutions  of  the  standard  wave  equation.  However,  the  method  is  not  very  efficient 
for  electrically  large  problems  because  one  has  to  solve  an  associated  huge  matrix  system. 
More  details  about  FEM  can  be  found  in  [5]. 

Different  from  MoM  and  FEM,  FDTD  is  a  time-domain  method.  The  currently  stan¬ 
dard  FDTD  method  for  electromagnetics  was  proposed  by  Yee  in  1966  [6].  In  Cartesian 
coordinates,  both  the  time  and  space  derivatives  in  Maxwell’s  equations  in  differential  form 
are  directly  replaced  by  second-order  central  difference  approximations.  Moreover,  Yee  clev¬ 
erly  staggers  the  positions  of  the  E  and  H  fields  to  naturally  model  the  physical  curls  of 
the  fields.  In  the  time  domain,  the  leap-frog  scheme  is  used  to  incrementally  march  the 


variation  of  the  E  and  H  fields.  Despite  the  simplicity  and  elegance  of  Yee’s  algorithm,  it 
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did  not  receive  much  interest  immediately  due  to  the  inability  to  model  an  ’’open”  problem 
and  the  lack  of  high  performance  computers.  As  the  computational  resources  became  more 
affordable,  FDTD  became  more  popular.  Since  the  late  80s,  the  number  of  publications 
related  to  the  FDTD  method  has  experienced  nearly  exponential  growth  in  [7].  Moreover, 
Berenger  presented  a  Perfectly  Matched  Layer  (PML)  [8]  -  [9] ,  which  successfully  truncated 
the  simulation  domain  and  increased  the  dynamic  range  of  FDTD  to  more  than  120  dB. 
Since  then,  FDTD  has  been  extensively  applied  to  almost  every  aspect  of  electromagnetics. 
As  a  numerical  method,  Yee’s  algorithm  is  explicit,  robust  and  easy  to  implement.  More 
importantly,  Yee’s  algorithm  is  able  to  model  more  complex  materials  than  any  other  nu¬ 
merical  methods  can  do.  However,  this  method  suffers  from  dispersion  errors  which  could 
contaminate  the  solutions  with  nonphysical  artifacts.  A  comprehensive  introduction  on  the 
Yee  algorithm  can  be  found  in  [10]  -  [14]. 

The  continuous  boost  of  the  information  industry  during  the  last  few  decades  has 
profoundly  changed  the  daily  life  of  everyone.  The  demands  for  better  information  prod¬ 
ucts  and  services,  such  as  high-definition  televisions,  the  future  wireless  LAN,  and  the 
third-generation  mobile  communications,  have  been  pushing  the  system  operating  frequency 
higher  due  to  the  following  reasons  [3].  First  of  all,  the  lower  end  of  the  electromagnetic 
spectrum  are  being  rapidly  depleted.  The  pursuit  of  higher  bit  rate,  which  requires  more 
bandwidth,  can  only  be  achieved  at  higher  frequencies.  Secondly,  the  system  miniaturiza¬ 
tion  requires  a  smaller  physical  size  of  the  antennas.  However,  the  gain  and  bandwidth  of 
an  antenna  are  proportional  to  its  electrical  size.  Therefore,  high  gain  and  bandwidth  are 
only  possible  at  a  higher  frequency.  As  a  consequence  of  the  general  industry  trend,  the 


operating  frequency  will  continuously  increase. 


5 


A  higher  operating  frequency  means  that  the  size  of  the  platform  that  carries  the 
devices,  such  as  a  mobile  vehicle,  a  commercial  or  military  aircraft,  becomes  electrically 
larger.  This,  however,  presents  challenges  to  the  numerical  methods  due  to  the  following: 

•  The  simulations  of  the  electrically  large  problems  require  to  store  vast,  sometimes 
prohibitive,  amount  of  variables  in  computer  memories.  For  example,  for  the  FDTD 
algorithms,  the  memory  needed  to  store  the  electric  and  magnetic  fields  is  proportional 
to  the  product  of  the  number  of  cells  in  each  direction  (Nx  x  Ny  x  Nz). 

•  Even  if  one  can  afford  the  demand  of  the  astonishing  amount  of  memory,  the  execution 
time  to  run  such  a  simulation  for  electrically  large  problems  is  extremely  long. 

•  Even  if  one  is  patient  enough  to  wait  until  the  completion  of  the  simulation,  chances 
are  the  results  are  totally  contaminated  by  the  numerical  errors,  especially  the  phase 
errors  due  to  the  dispersion  . 

The  currently  standard  EDTD  algorithm  (Yee’s  scheme)  is  only  second-order  accurate 
in  both  the  space  and  time  domains,  and  suffers  from  severe  nonphysical  dispersions.  Based 
on  many  analytical  analysis  and  numerical  experiments  [1],  it  has  been  found  that  a  cell  size 
of  A/10  —  A/20  provides  acceptable  accuracy  for  most  problems  with  moderate  electrical 
size.  However,  the  phase  errors  generated  by  dispersion  accumulate  as  the  traveling  distance 
of  the  wave  and  the  simulation  time  increase.  Consequently,  for  electrically  large  problems, 
this  rule  of  thumb  may  not  be  valid  since  the  accumulation  of  the  phase  errors  due  to  the 
dispersion  significantly  deteriorate  the  accuracy.  A  simple  way  to  reduce  the  dispersion 
errors  is  to  use  a  finer  mesh  resolution.  This  yields  an  even  larger  simulation  domain 


which  aggravates  the  already  heavy  burden  of  the  computer  memory  and  CPU.  In  other 
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words,  the  available  computer  resources  limit  the  electrical  size  of  the  problems  that  can  be 
handled.  Hence,  the  mesh  refinement  is  not  an  efficient  solution  and  sometimes  not  even 
an  alternative. 

In  reality,  the  dispersion  is  one  of  the  inherent  characteristics  related  to  the  FDTD 
methods  which  cannot  be  completely  eliminated.  In  order  to  minimize  the  dispersion, 
many  alternative  approaches  were  proposed.  These  methods,  in  the  broadest  sense,  can 
be  classified  into  two  categories:  the  Standard  FDTD  (SFD)  and  the  Nonstandard  FDTD 
(NSFD).  The  standard  FDTD  schemes  are  derived  based  upon  the  Taylor  series  expansion. 
The  coefficients  of  the  stencils  are  chosen  to  minimize  the  truncation  error  and  maximize 
the  order  of  accuracy.  As  a  comparison,  the  nonstandard  FDTD  schemes  consider  the 
minimization  of  dispersion  as  the  ultimate  goal.  Instead  of  using  the  Taylor  series  expansion, 
the  coefficients  of  the  NSFD  stencils  are  derived  directly  from  the  dispersion  relations. 

Yee’s  algorithm  is  a  typical  standard  FDTD  scheme  with  second-order  accuracy  in  both 
the  time  and  space  domains.  Conventionally,  an  SFD  method  with  Mth-  and  Nth-order 
accuracy  in  the  time  and  space  domains,  respectively,  is  referred  to  as  the  FDTD(M,N) 
scheme.  For  example,  Yee’s  algorithm  can  be  also  denoted  as  the  FDTD(2,2)  scheme.  It 
is  a  natural  consideration  to  attempt  to  reduce  the  dispersion  by  increasing  the  stencil’s 
order  of  accuracy.  In  1989,  Fang  presented  a  scheme  which  is  second-order  accurate  in 
the  time  domain  but  fourth-order  accurate  in  the  space  domain  [15].  Fang’s  algorithm  is 
also  referred  to  as  the  FDTD(2,4)  scheme.  Since  then,  higher-order  schemes  quickly  drew 
the  attention  of  the  FDTD  researchers  and  became  one  of  the  most  promising  dispersion- 
reduction  methods.  Early  exploration  of  Fang’s  algorithm  and  other  high-order  FDTD 
schemes  includes  [16]-  [25].  A  modified  FDTD(2,4)  scheme  based  on  the  integral  form  of 
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Maxwell’s  equations  was  introduced  in  [26],  [27].  A  second-order  “isotropic”  scheme  using 
complex  multi-dimensional  stencils  was  proposed  in  [28].  Contributions  in  the  time-domain 
discretization  include:  implicit  time- integration  schemes  [29],  [30]  and  higher-order  time- 
integration  schemes  [31]  -  [32] .  Some  intuitive  reviews  and  computational  efficiency  analysis 
of  higher-order  FDTD  can  be  found  in  [33]  -  [36]. 

Another  important  type  of  method  is  the  nonstandard  FDTD  (NSFD).  The  stencils  of 
the  NSFD  may  be  similar  to  the  stencils  of  the  SFD.  However,  the  NSFD  coefficients  are 
properly  adjusted  from  their  SFD  counterparts  for  further  dispersion  reduction.  Typically, 
the  values  of  the  NSFD  coefficients  are  determined  by  the  Fourier  analysis  based  on  the  dis¬ 
persion  relation.  The  concept  of  the  NSFD  was  first  proposed  by  Mickens  [37].  Applications 
of  the  NSFD  in  various  mathematical  and  physical  areas  are  summarized  in  [38],  [39].  Cole 
did  the  pioneering  work  in  applying  the  NSFD  to  solve  the  Maxwell’s  equations  [40],  [41].  To 
minimize  the  anisotropy  of  the  NSFD  for  multi-dimensions.  Cole  proposed  a  complex  stencil 
which  is  similar  to  the  one  used  in  [28].  Afterwards,  the  NSFD  scheme  is  modified  to  be  able 
to  model  the  lossy  and  absorbing  materials  [42],  [43].  A  full  dispersion  and  stability  analysis 
for  the  3D  Cole’s  scheme  can  be  found  in  [44] ,  [45] .  Furthermore,  Cole’s  scheme  is  expanded 
to  even  more  complex  stencils  and  curvilinear  coordinate  systems  in  [46]-  [52].  In  order 
to  reduce  the  average  dispersion  error  without  using  the  complex  stencils,  [53]-  [55]  uses 
dispersion  centering  technique  originated  in  [1].  Then,  a  scheme  that  can  control  the  disper¬ 
sion  level  in  a  certain  range  of  propagation  angles  was  presented  by  introducing  a  filtering 
process  [56],  [57].  Although  proposed  independently  by  different  research  groups,  [58]-  [65] 
indicated  that  the  dispersion  error  can  be  drastically  reduced  by  using  degenerated  fourth- 


order  stencils.  That  is,  the  coefficients  of  the  fourth-order  stencil  are  deliberately  adjusted 


to  further  reduce  the  dispersion.  The  cost  is  that  the  accuracy  of  the  fourth-order  stencil  is 
degenerated  to  second-order.  Comparisons  of  the  dispersion  properties  of  Cole’s  and  other 
low-dispersion  algorithms  can  be  found  in  [66],  [67]. 

In  order  to  further  decrease  the  dispersion  error,  the  previously  stated  standard  and 
nonstandard  FDTD  schemes  both  use  extended  stencils.  The  extended  stencils  represent 
those  finite  difference  operators  that  involve  more  than  two  points.  For  example.  Fang’s 
FDTD(2,4)  scheme  is  considered  to  be  a  simple  extended  stencil  with  four  points.  Typ¬ 
ically,  if  more  points  are  included  into  the  stencil,  more  coefficients  (degrees  of  freedom) 
can  be  used  to  mitigate  the  dispersion.  The  coefficients  of  the  extended  stencils  can  be 
derived  using  either  the  Taylor  series  expansion  (SFD)  or  dispersion  analysis  (NSFD).  In  a 
homogeneous  region,  the  schemes  using  the  extended  stencils  work  perfectly  well.  However, 
significant  nonphysical  errors  will  be  generated  at  the  material  interface.  These  errors  could 
be  attributed  to  the  field  discontinuities  along  each  side  of  the  material  interface.  Many 
efforts  have  been  devoted  to  restore  the  accuracy  of  the  extended  stencils  across  the  ma¬ 
terial  interface  [68]-  [71].  Moreover,  a  subgridding  technique  to  hybridize  the  FDTD(2,2) 
and  the  FDTD(2,4)  [72]  -  [76]  was  introduced.  The  FDTD(2,4)  scheme  is  only  used  in 
the  homogeneous  regions.  However,  along  the  material  interfaces  and  the  domain  bound¬ 
aries,  the  FDTD(2,2)  scheme  is  applied.  Along  the  interface  between  the  FDTD(2,2)  and 
FDTD(2,4)  domains,  the  fields  are  carefully  coupled.  This  method  does  not  require  any 
special  treatment  for  the  material  interface. 

All  of  the  conventional  FDTD(2,2)  truncation  techniques,  such  as  the  PML,  can  be 
used  without  modifications.  However,  the  hybrid  scheme  is  still  second-order  accurate 


and  the  programming  is  cumbersome.  Alternative  approaches  include  the  one-sided  stem 
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cils  [77],  [78].  The  derivatives  near  the  material  interface  are  approximated  using  points  only 
along  one  side  of  the  interface.  However,  as  reported  in  [77] ,  the  accuracy  of  the  one-sided 
stencil  degenerates  quickly  when  the  discontinuity  is  aggravated.  Originated  in  [79],  [80],  a 
derivative  matching  method  was  proposed  [81].  This  method  introduces  fictitious  points  at 
both  sides  of  the  material  interface  and  physically  enforces  the  boundary  conditions.  Using 
this  technique,  consistent  stencils  can  be  applied  throughout  the  entire  domain.  In  [81],  it 
has  been  shown  that  the  accuracy  of  the  schemes  (up  to  eighth-order)  are  fully  restored. 
Another  method  [82]-  [84]  is  based  on  the  Fourier  analysis  of  the  numerical  reflection  co¬ 
efficient  which  is  very  similar  as  the  derivation  of  the  NSFD  schemes.  Although  originally 
presented  for  Yee’s  algorithm,  this  method  seems  very  promising  for  the  extended  sten¬ 
cils.  Useful  reviews  on  the  state-of-art  progresses  in  the  extended  stencil  and  the  material 
interface  treatments  can  be  found  in  [85],  [86]. 

It  has  been  indicated  that  the  dispersion  can  be  reduced  significantly  by  using  the 
NSFD  methods.  However,  the  following  technical  issues  still  need  to  be  considered. 

•  For  multi-dimensional  problems,  the  NSFD  schemes  still  suffer  from  serious  nonphys¬ 
ical  anisotropic  errors.  Methods  to  decrease  the  anisotropy  without  significantly  in¬ 
creasing  the  computational  burden  are  still  highly  desired. 

•  The  current  available  approaches  to  determine  the  values  of  the  stencil  coefficients 
become  very  cumbersome  for  more  complex  schemes.  Simple  approaches  to  solve  the 
optimal  coefficients  are  still  needed. 

•  Some  principal  concerns  for  the  NSFD  algorithms,  such  as  the  Gauss’s  Law,  have  not 


been  examined  yet. 
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•  The  available  techniques  treating  the  material  interface  are  far  from  practical.  Ironi¬ 
cally,  it  is  usually  much  easier  to  design  a  very  complex  stencil  than  the  corresponding 
interface  conditions.  Proper  material  interface  conditions  are  still  the  greatest  chal¬ 
lenge  in  this  area. 

With  the  increase  of  the  operating  frequency,  electromagnetic  compatibility  (EMC) 
and  electromagnetic  interference  (EMI)  issues  become  incrementally  more  important.  One 
of  the  most  critical  EMI  concern  is  that  the  aircraft  navigation/communication  antennas 
could  be  interfered  by  the  radiation  emission  generated  by  on-board  personal  electronic 
devices  (PEDs),  such  as  cellphones,  laptops,  etc.  That  is  the  reason  for  switching  off  all 
the  personal  electronic  devices  during  the  take-off  and  landing,  if  not  during  the  entire 
flight.  This  EMI  issue  has  be  studied  by  characterizing  the  mutual  coupling  between  an 
on-board  PED  antenna  and  an  aircraft  antenna  mounded  on  the  exterior  of  a  simplified 
Boeing  757  fuselage  [72].  However,  the  fuselage  cavity  in  the  previous  research  was  empty. 
In  practice,  the  fuselage  is  always  filled  with  various  metallic  frames,  decoration  materials 
and  passengers.  Among  these  materials,  the  human  body  could  be  the  most  important  due 
to  its  highly  dielectric  and  lossy  property.  The  impact  of  the  present  passengers  on  the  EMI 
characteristics  still  remains  a  question. 

In  this  project,  the  major  objective  is  to  resolve  or  at  least  mitigate  the  aforementioned 
technical  problems  for  the  nonstandard  EDTD  methods.  The  passenger  effects  on  the  PED 
mutual  coupling  will  also  be  examined.  The  structural  organization  of  this  report  is  now 
outlined.  In  Chapter  2,  the  fundamental  principles  of  the  electromagnetics  will  be  first 
reviewed.  Then  the  basic  theory  of  the  standard  EDTD,  including  the  EDTD(2,2)  (Yee’s) 
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and  other  higher-order  schemes,  will  be  presented.  A  detailed  dispersion  analysis,  which 
forms  the  basis  of  the  nonstandard  FDTD,  on  various  standard  FDTD  algorithms  will  also 
be  discussed. 

Chapter  3  is  focused  on  the  nonstandard  FDTD.  We  will  go  through  the  fundamental  of 
the  one-dimensional  NSFD  algorithm.  Dispersion  analysis  and  numerical  simulations  will  be 
presented  to  reveal  the  strength  and  weakness  of  the  ID  NSFD.  An  NSFD(2,4)  scheme  will 
also  be  introduced  by  incorporating  the  concept  of  NSFD  into  the  fourth-order  stencil.  The 
NSFD  algorithms  will  then  be  extended  to  multi-dimensional  cases.  To  reduce  the  average 
dispersion  error,  an  improved  NSFD  algorithm  will  be  presented.  A  simple  least  squares 
method  will  be  used  to  determine  the  optimal  stencil  coefficients.  Afterwards,  the  concept 
of  the  NSFD  will  be  generalized.  Dispersion  analysis  will  be  presented  to  demonstrate  the 
performance  of  the  generalized  NSFD.  Then,  the  Gauss’s  Laws  in  both  the  differential  and 
integral  forms  will  be  examined  for  the  NSFD  algorithms.  Finally,  a  material  interface 
condition  based  on  the  concept  of  the  NSFD  will  be  proposed  to  decrease  the  errors  caused 
by  the  material  interface. 

In  Chapter  4,  the  passenger’s  effect  on  the  FED  mutual  coupling  will  be  investigated. 
The  mutual  coupling  between  the  antennas  is  described  by  the  5-parameters.  A  simplified 
passenger  model  was  designed  and  built  to  fit  into  the  fuselage  model.  The  5-parameters 
will  be  simulated  using  the  FDTD  algorithm.  To  justify  the  simulations,  same  parameters 
are  measured  using  the  vector  network  analyzer  in  the  electromagnetic  anechoic  chamber 
of  Arizona  State  University.  Finally,  the  5-parameters  with  and  without  passengers  are 


compared  to  demonstrate  the  passenger’s  effect. 


CHAPTER  2 


THE  STANDARD  EINITE-DIEEERENCE  METHODS 

As  indicated  in  the  precious  chapter,  the  Standard  Einite-Difference  Time-Domain 
Method  (SED)  is  an  important  type  of  EDTD  schemes  that  is  derived  directly  using  the 
Taylor  series  expansion.  The  simplest  SED  method,  Yee’s  algorithm,  is  the  currently 
most  popular  EDTD  scheme  that  has  been  extensively  applied  to  many  electromagnetic 
boundary- value  problems,  for  example,  antenna  modeling,  radiation/scattering,  microwave 
circuits  design,  and  EMC/EMI  analysis,  etc.  However,  Yee’s  algorithm  is  only  second-order 
accurate  in  both  the  time  and  space  domain.  Hence  it  suffers  from  severe  dispersion  error 
especially  for  the  electrically  larger  problems.  In  order  to  reduce  the  dispersion  error,  one  of 
the  most  promising  methods  is  to  use  the  higher-order  SED  method  (HOSED).  Usually,  a 
central-difference  space  stencil  using  N  points  could  achieve  an  order  of  accuracy  up  to  A^-th 
order.  However,  the  accuracy  of  the  entire  scheme  also  depends  on  the  time-domain  dis¬ 
cretization  and  the  type  of  the  problems.  Moreover,  due  to  the  discontinuities  of  the  fields, 
significant  errors  can  be  generated  along  the  material  interface  which  will  deteriorate  the  ac¬ 
curacy  of  the  entire  HOSED  method  to  be  similar  as  that  of  Yee’s  algorithm  [85],  [86].  Since 
increasing  the  number  of  points  doses  not  achieve  the  desirable  accuracy,  many  advances 
on  the  HOSED  only  consider  Eang’s  EDTD (2, 4)  scheme,  which  is  second-order  accurate  in 
the  time  domain  but  fourth-order  accurate  in  the  space  domain. 

In  this  chapter,  the  fundamental  principles  of  any  numerical  method,  i.e.  Maxwell’s 
equations,  will  be  first  reviewed.  Then  the  theory  of  Yee’s  algorithm  will  be  introduced.  A 
general  procedure  to  derive  the  SED  stencils  with  high  order  of  accuracy  will  follow.  The 


stability  conditions  for  the  SED  methods  will  be  discussed.  Einally,  a  detailed  dispersion 
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analysis  for  Yee’s  and  Fang’s  algorithms,  which  are  also  the  base  of  the  nonstandard  FDTD 
methods,  will  be  presented. 

2.1.  Fundamental  Principals 

Time-varying  electromagnetic  phenomena  are  initial  boundary-value  problems.  In  a 
homogeneous  region,  these  problems  are  fully  described  by  Maxwell’s  equations.  Along  the 
material  discontinuities,  proper  boundary  conditions  must  be  enforced.  Maxwell’s  equations 
can  be  written  either  in  differential  or  in  integral  form. 

The  differential  form  of  Maxwell’s  equations  is  most  widely  used  since  it  governs  the 
relations  between  field  vectors,  current  densities,  and  charge  densities  at  any  point  in  space 
at  any  time.  For  a  homogeneous  region.  Maxwell’s  equations  in  differential  form  can  be 
written  as 


VxE  =  — 
ot 

(2.1) 

<9D 

V  X  H  —  — - — h  Jc  +  Jj 
ot 

(2.2) 

V  •  D  =  ge 

(2.3) 

V  •  B  =  q^n 

(2.4) 

where  E  is  the  electric  field  intensity  in  (volts/meter),  H  is  the  magnetic  field  intensity 
in  (amperes/meter),  D  is  the  electric  flux  density  in  (coulombs/square  meter),  B  is  the 
magnetic  flux  density  in  (webers/square  meter),  Jc  and  Jj  are  the  conduction  and  impressed 
electric  current  densities,  respectively,  in  (amperes /square  meter).  Me  and  Mj  are  the 
conduction  and  impressed  magnetic  current  densities,  respectively,  in  (volts/square  meter), 
qe  is  the  electric  charge  density  in  (coulombs/cubic  meter),  and  qm  is  the  magnetic  charge 
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density  in  (webers/cubic  meter).  Individually,  (2.1)-(2.4)  are  also  referred  to  as  Faraday’s 
Law,  Ampere’s  Law,  Gauss’s  Law  for  electric  field,  and  Gauss’s  Law  for  magnetic  field, 
after  the  names  of  their  discoverers. 

The  integral  form  of  Maxwell’s  equations  describes  relations  of  the  field  vectors,  current 
densities,  and  charge  densities  over  an  extended  region  of  space.  The  integral  form  of 
Maxwell’s  equations  can  be  written  as 


Edl  = 


c 


Hdl 


c 


(2.5) 

(2.6) 

II 

Q 

(2.7) 

^  B  ■  ds  =  Qm 

(2.8) 

where  Qe  and  Qm  are  the  total  electric  and  magnetic  charges,  respectively. 

In  materials  with  constitutive  parameters  that  are  independent  of  time,  the  electric  and 
magnetic  flux  densities  (D  and  B)  are  related  to  the  electric  and  magnetic  field  intensities 
E  and  H  using  the  constitutive  relations 


D  =  eE 

(2.9) 

B  =  /iH 

(2.10) 

where  e  is  the  electrical  permittivity  of  the  medium  in  (Farad/meter)  and  fx  is  the  magnetic 
permeability  of  the  medium  in  (Henries/meter).  For  free  space,  e  and  jj,  are  constants  which 
can  be  written  as 


eo  =  8.854  x  10  {Farads /meter) 


(2.11) 
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/io  =  47r  X  10  {H envies / meter)  (2-12) 

Moreover,  the  conduction  electric  current  density  Jc  and  the  conduction  magnetic 
current  density  Me  are  related  to  the  electric  and  magnetic  field  intensities  E  and  H  using 
the  Ohm’s  laws 


J  =  cjE 

(2.13) 

M  =  u*H 

(2.14) 

where  a  is  the  electrical  conductivity  in  (Siemens/meter)  and  a*  is  the  magnetic  conduc¬ 
tivity  in  (Ohms/meter). 

In  order  to  completely  describe  the  electromagnetic  phenomena,  appropriate  boundary 
conditions  must  be  enforced  along  the  interfaces  where  the  media  exhibit  discontinuities. 
Assuming  there  are  no  sources  present,  along  a  interface  formed  by  two  media  with  electric 
properties  (ei,/ri)  and  (€2,1^2),  the  boundary  conditions  can  be  expressed  as 


h  X  (El  -  E2) 

=  0 

(2.15) 

h  X  (Hi  -Hs) 

1  =  0 

(2.16) 

h  ■  (Di  -  D2) 

=  0 

(2.17) 

h  •  (Bi  -  B2) 

=  0 

(2.18) 

where  h  is  the  unit  vector  normal  to  the  interface,  pointing  from  medium  2  into  medium 
1.  A  very  important  special  case  is  that  medium  2  is  a  perfect  electric  conductor  (PEC). 
Since  the  fields  inside  the  PEC  are  zero,  (2.15)-(2.16)  reduce  to 


h  X  El  =  0 


(2.19) 


16 


n  X  Hi  =  0  (2.20) 

That  is,  along  the  PEC  interface,  the  tangential  components  of  the  electric  and  magnetic 
fields  are  zero. 

The  two  curl  equations  of  (2.1)-(2.2)  form  the  basis  of  the  FDTD  algorithm  for  solving 
3D  electromagnetic  problems.  Usually,  the  two  Gauss’s  Laws  of  (2.3)-(2.4)  will  not  be 
explicitly  implemented  in  FDTD.  However,  to  avoid  the  risk  of  introducing  artificial  charges 
into  the  domain,  it  is  very  important  to  check  whether  or  not  the  two  Gauss’s  laws  are 
implicitly  satisfied  by  any  FDTD  algorithm. 


2.2.  The  Standard  FDTD(2,2)  (Yee’s)  Algorithm 


As  the  simplest  case  of  the  SFD  schemes,  the  FDTD(2,2)  (Yee’s)  algorithm  is  second- 
order  accurate  in  both  the  time  and  space  domains.  In  other  words,  the  time  and  space 
derivatives  in  Maxwell’s  equations  are  replaced  by  the  second-order  central  difference  ap¬ 
proximations.  Letting  /(^,  t)  represent  an  arbitrary  function  of  space  and  time,  this  can  be 
expressed  as 


dt 


/(C +  Ag/2,t)-/(e- AC/2,  t) 

AC 

/(C,  t  +  At/2) -/(C,t  + At/2) 
At 


+  0(AC'), 

+  O  (At^) 


i  =  x,y,z 


(2.21) 

(2.22) 


Gonventionally,  the  finite  difference  approximations  (2.21)-(2.22)  are  also  referred  to 
as  the  space  and  time  stencils^  respectively.  For  example,  the  space  stencil  along  the  x-axis 
is  illustrated  in  Fig.  2.1.  As  shown,  in  order  to  evaluate  the  partial  derivative  df  /dx  at  the 
center  point  x,  one  needs  to  calculate  the  difference  between  the  two  neighboring  values  at 
X  —  Ax/2  and  x  -|-  Ax/2;  then  divide  it  by  the  distance  between  them  (Ax).  The  constants 
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Fig.  2.1.  Second-order  central  difference  stencil  along  the  x  direction. 


above  each  point  are  the  weighting  coefficients  which  can  be  derived  from  Taylor  series 
expansion. 

Before  the  introduction  of  Yee’s  algorithm,  the  notations  that  follows  will  be  introduced 
and  used  repeatedly.  In  Cartesian  coordinates,  the  space  steps  with  respect  to  x,  y,  and  z 
directions  and  the  time  increment  are  denoted  as  Ax,  Ay,  Az,  and  At,  respectively.  The 
space  position  at  {iAx,jAy,kAz)  is  represented  by  the  space  indices  {i,j,k).  The  time 
instant  nAt  is  represented  by  the  time  index  n.  Following  these  notations,  an  arbitrary  grid 
function  at  time  t  =  nAt  can  be  expressed  as  f\^jk- 

To  demonstrate  Yee’s  algorithm,  let  us  assume  that  a  plane  wave  is  traveling  in  the  z 
direction  in  an  unbounded,  homogeneous  medium  (e, /r,  ct,  a*).  Without  loss  of  generality, 
we  also  assume  that  the  plane  wave  has  an  component  and  a  magnetic  field  component 
Hy.  as  illustrated  in  Fig.  2.2. 

In  Cartesian  coordinates,  the  plane  wave  can  be  modeled  using  the  one-dimensional 
Maxwell’s  equations  which  can  be  written  as 


dE^ 

dt 


e 


dz 


+  <7  Ex 


(2.23) 
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=  Ca{k)  -  C,{k){Hy\l^y2  -  Hy\l_,/^)  (2.25) 

Hy\lX\,2  =  Da{k  +  l/2)Hy\l^^j^-Dy{k  +  l/2){EXuX\'^-ESk""'^^  (2.26) 


where 


Ca{k) 

Cb{k) 

Da{k) 

Db{k) 


1  -  (cjfcAt)/(2efc) 
1  +  (cjfcAt)/(2efc) 

At/{ekAz) 

1  +  (cjfcAt)/(2efc) 
1  -  (a^At)/(2^fc) 
1  +  (c7*At)/(2^fc) 
At/inkAz) 

1  +  ((T*At)/(2^fc) 


(2.27) 

(2.28) 

(2.29) 

(2.30) 
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Fig.  2.3.  Time-space  chart  for  the  one-dimensional  Yee  algorithm. 

are  the  update  coefficients.  Notice  that  e  and  a  take  the  values  at  the  positions  where  the 
E  fields  are  located  and  /x  and  a*  take  the  values  where  the  H  fields  are  located. 

Yee’s  algorithm  can  be  also  explained  using  a  time-space  chart  illustrated  in  Fig.  2.3 
[13].  As  shown,  to  update  the  field  at  (A:, n -|- 1/2),  one  needs  the  E^  field  at  (A:,n  — 1/2) 
and  the  Hy  fields  at  {k  —  1/2,  n)  and  (/c  -|-  1/2,  n).  After  that,  the  updated  Ex  held  at 
(/c,  n  1/2),  together  with  the  Ex  held  at  {k-\-l^n-\- 1/2)  and  the  Hy  field  at  {k  -|-  1/2,  n), 
will  be  used  to  update  the  Hy  field  at  {k  +  1/2,  n  -|-  1).  This  progress  continues  until  the 
time  advance  is  terminated. 

In  two  dimensions,  electromagnetic  boundary- value  problems  can  be  classified  by  two 
modes;  i.e.,  TE^  and  TM^  modes.  In  Cartesian  coordinates.  Maxwell’s  equations  for  TE^ 


20 


(2.31) 

(2.32) 

(2.33) 

(2.34) 

(2.35) 

(2.36) 


Staggering  the  electric  and  magnetic  fields  in  both  the  x  and  y  directions,  the  two- 
dimensional  Yee’s  grids  for  TE^  and  TM^  modes  are  demonstrated  in  Figs.  2.4  and  2.5, 
respectively. 

In  the  time  domain,  the  leap-frog  scheme  is  used  to  discretize  the  time  derivatives. 
Following  a  similar  procedure  as  for  one-dimensional  problems,  the  two-dimensional  Yee’s 
algorithms  for  and  TE^  and  TM^  modes  take  the  form  of  [6],  [1] 
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TE^  mode: 

=  cji + 1/2.  j) 

+  +  l/2,i)(-f^2|j+i/2j+l/2  “  -^z|i+i/2J-l/2)  (2.37) 

=  C^i,  j  + 1/2) 

—  Cbxii,  j  +  ^/‘^)iHz\^_^_l/2,j+l/2  ~  ^z\i-l/2,j+l/2)  (2.38) 

^z\i^l/2,j+l/2  ~  Da{i  +  1/2,  j  +  1/2)  Hz\^_^_i/2,j+l/2 

+  Di,{i  +  1/2,  j  +  l/2){E^\’;+;ll+i  -  EA"+;ll) 

-  Dt.b  +  l/2,j  +  l/2)(.E,\’;+;;i‘_^,,,- E,\'l*l{%)  (2.39) 

TM^  mode: 

Da{i  +  1/2,  j)  77a;|j+i/2j 

+  1/2,  -  E,isyi._v2)  (2.40) 

Da{i,j  +  1/2)  .f^i/ljj+1/2 

Dbx{i,j  +  l/2)(^z|i+l/2j+i/2  “  -®2|j_i/2J+l/2)  (2-41) 

Ca{i  +  1/2,  j  +  1/2) 

Cby{i  +  1/2,  j  +  l/2)(77a;|j+i/2,j+l  “  -^oi|i+i/2j) 
a,(z  +  l/2,j  +  l/2)(F,|/^,^^.^,/2-  ^yi::,+i/2)  (2-42) 


TJ  \Ti-\-l  _ 


TJ  ITI+I  _ 


p  |?i.+l/2  _ 

^^\i+l/2,j+l/2  - 
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where 


Ca{i,j) 

Cbdhj) 

Da{i,j) 

Dbdhj) 


1  -  {aijAt)/{2eij) 

1  +  /{2eij) 

At/{eijAC)  ^ 

1  +  (fTjjAt)/(2ejj’ 

1  -  «^.At)/(2^,,,) 

1  +  {al^At)/{2^^iJ) 
At/{iMijA() 
l  +  {al^At)/{2^i,j)’ 


(2.43) 

(2.44) 

(2.45) 

(2.46) 


In  the  three-dimensional  Cartesian  coordinates,  Maxwell’s  equations  of  (2.1)  and  (2.2) 
reduce  to  a  system  of  six  coupled  scalar  equations  which  can  be  expressed  as  [2] 


dHr 


dH,, 


OH,  ^ 
dt  i^y  dy 
dE^  _  1 
dt  e 

dEy  _  1 


dt 

dE, 

dt 


dEy 

dE, 

dz 

dy 

dE, 

dE, 

dx 

dz 

dE, 

dEy 

dy 

dx 

dH, 

dHy 

dy 

dz 

dH, 

dH, 

dz 

dx 

dHy 

dH, 

dx 

dy 

-a*H, 


-  a*H, 


-  a*H, 

—  aEr 


—  aE, 


—  aE. 


(2.47) 

(2.48) 

(2.49) 

(2.50) 

(2.51) 

(2.52) 


To  solve  three-dimensional  problems,  the  space  domain  is  discretized  using  the  Yee’s  unit 
cell,  as  illustrated  in  Fig.  2.6.  As  shown,  all  the  field  components  are  located  at  different 
positions  on  the  surface  of  the  Yee  cell.  The  E  and  H  fields  are  staggered  at  a  distance 
of  half  space-step  which  allows  the  explicit  evaluation  of  the  curls  at  each  surface  plane. 
Different  constitutive  parameters  are  locally  assigned  to  each  field  component  to  model 


different  materials.  Using  the  leap-frog  scheme  in  the  time-domain,  the  three-dimensional 
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Fig.  2.6.  The  three-dimensional  Yee’s  unit  cell  [1]. 


Yee  update  equations  are  given  by  [6] ,  [1] 


+  -|-  1/2,  j,  fc)(i7^|^_,_j^y2j+l/2,fc  “  -^2:|j-|-l/2j-l/2,fc) 

—  Cbz{i+  1/2,4,  A:)(i?j;|j_,_j^y2j,fc+l/2  “  -^yli+i/2j,fc-i/2) 
^y\i,j+i/2,k  ~  +  1/^’ -®?/liJ+l/2,fc 

+  Cbziijj  +  ^/‘^,k){Hx\  i,j+l/2,k+l/2  ^^\i,j+l/2,k-l/2) 
—  Cbxihj  +  1/2,  k){Hz\i_^_i/2,j+l/2,k  ~  ^z\i-l/2,j+l/2,k) 

^•eSS/2  =  c'„{i.;.*,+i/2)E,i"-;S/2 

+  Cbx{i,j,k  +  l/2)(-f^2;|j_,_]^/2j,fc+l/2  “  ^y\i-l/2,j,k+l/2') 
—  Cby{i,j,k  +  l/2)(i^x|ij+i/2,fc+l/2  ~  -l/2,k+l/2) 


Cbzii  +  - 


(2.53) 


■  ^z\i-l/2,j+l/2,k)  (2-54) 


(2.55) 
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where 


TT  \n-\-l 


TT  I  n.+  l 

^y\i+l/2,j,k^l/2 


=  Da{i,  j 1/2,  k 1/2)  Hx\ij^i/2,k+l/2 
+  Djjzihj  +  1/2,  ^  +  l/2)(-E'y|j^^_,_(/2,fc+i 
—  Df,y{i,j  +  1/2,  k  +  1/‘2){Ez\-j_^{^i._^-^^2 
Da{i  +  l/2,j,k+l/2) 


-  E, 


-  E. 


|n+l/2  \ 

'yHJ+l/2,k) 

|n+l/2  \ 

^l^,_7,/c+l/2■' 


I  ^.ir+ 


n 

i+l/2J,k-\-l/2 

+  DUi  +  1/2,  j,  +  1/2) ( 


JT  I 


+ 


I  ii'-ri-j  z 

H+l,j,k+l/2  ' 

Dbzii  +  l/‘2,j,k+  l/2)(-E'a;|^_l_^y2j,A:+l  ' 

+  1/2,  J  +  1/2,  k)  ^^2:|i+l/2j  +  l/2,fc 
Dby{i  +  1/2,  J  +  1/2, 

Dbx{i  +  l/2,i  +  1/2,  ^)(-E'ylj+i  j+i/2,fc 


|n+l/2  \ 

I  h+1/2 


^U+l/2j,/c-' 


P  |n+l/2  X 

^U+l/2,j,fc-' 
|n+l/2  \ 


- 


Ca{i,j,k)  = 
Cb^{hj,k)  = 
Da{i,j,k)  = 
Dbi{i,j,k)  = 


1  (^j,j,fe^^)/(2ei  j,fc) 

1  “1“  (o'jjjfcAt)/(2ej 

At/(e^,j-fcAO 

1  +  (o-jj,fcAt)/(2eij-fc’ 
1  ~  ^^i,j,k^^) ! ^‘^k'i,j,k) 
1  +  /  ^‘^k'i,j,k) 

At/(/Ujj,fcA^) 


i  =  x,y,z 


,  i  =  x,y,z 


(2.56) 


(2.57) 


(2.58) 


(2.59) 

(2.60) 
(2.61) 
(2.62) 


1  (^i,j,k‘^^)/ ^‘^k'i,j,k) 

Despite  its  second-order  accuracy,  Yee’s  algorithm  is  the  most  popular  FDTD  method 
for  computational  electromagnetics.  For  years,  it  has  been  successfully  applied  to  solve  nu¬ 
merous  practical  electromagnetic  problems.  There  is  a  large  collection  of  models  simulated 
by  Yee’s  grid,  from  a  small  IC  package  to  a  large  fuselage.  When  a  dispersion-reduction 
method  is  designed,  it  is  usually  desired  that  the  new  algorithm  is  able  to  take  advantage 
of  the  available  geometry  library.  In  other  words,  the  new  dispersion-reduction  scheme 
must  also  use  Yee’s  grid.  One  of  the  most  promising  approaches  is  to  use  the  standard 


26 


higher-order  schemes. 

2.3.  The  Standard  Higher-Order  Schemes 

The  standard  higher-order  schemes  denote  those  methods  which  use  the  finite-difference 
approximations  with  accuracy  higher  than  second-order.  Higher-order  stencils  can  be  used 
to  discretize  the  derivatives  in  both  the  time-  and  the  space-domain.  For  clarity,  we  only 
consider  the  space-domain  higher-order  stencils.  In  the  time  domain,  the  second-order 
stencil  is  still  used.  Hence,  the  standard  higher-order  methods  in  this  report  can  be  referred 
to  as  the  FDTD(2,A^)  schemes.  Notice  that,  for  the  central  difference  stencils,  N  is  always 
an  even  integer.  To  be  compatible  to  the  the  available  geometry  file,  Yee’s  grids  are  used 
for  the  space  discretization. 

Before  the  discussion  of  the  characteristics  of  the  higher-order  stencils,  the  method  to 
derive  these  stencils  will  be  introduced.  For  example,  to  approximate  the  derivative  at  the 
center  point,  the  central  difference  approximation  involving  four  neighboring  points  can  be 
expressed  as 

Cif  (e  -  ^)  +  C2f  U-^)+  Csf  U  +  ^)+  c^f  (e  + 

/(O  ^ ^ ^ ^ ^ ^ - ^(2.63) 

C  =  x,y,z 


where  the  four  coefficients  C'i,C'2,C'3,  and  are  undetermined.  Expanding  the  function 
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values  at  the  four  neighboring  points  yields 


/U-^) 


fii-^ 


/le+f) 


/  (0  -  (0  +  ^/"  (f)  -  («)  +  (0 


/<"'  («) 


81A^5 
1280 

/  ({)  -  (0  +  ^/"  ({)  -  ^/"'  (9  +  ({) 

3840 


/(5)  (^)  +  ... 

AC 


(2.64) 


(2.65) 


/  (0  +  f^/'  (9  +  ^/"  (9  +  (9  +  (9 


A^3 


+ 


A^5 

3840 


(C)  + 


(2.66) 


/  (9  +  ^/'  (9  +  ^/"  (9  +  ^/"'  (9  +  (9 


+ 


81A^5 

1280 


(C) 


(2.67) 


Substituting  (2.64)-(2.67)  into  (2.63)  and  collecting  terms  leads  to 

/'(C)  (Ci  +  C2  +  C3  +  C4)^/(C) 

+  (_3Ci-C2  +  C3  +  3C4)2/'(C) 

+  (OCi  +  (72  +  Cs  +  9(74) -^/"  (^) 

+  (_27Ci-C2  +  C'3  +  27C4)^/'"(C) 

+  (81(7i  +  (72  +  (73  +  81(74)-^^/^^^  (^) 

+  (-243Ci-C2  +  C3  +  243C4)^/(^)(C) 

+  •  •  •  (2.68) 


In  order  to  approximate  the  derivative,  the  coefficient  of  the  second  term  in  (2.68)  must 


be  unity.  Furthermore,  to  achieve  the  accuracy  as  high  as  possible,  the  coefficients  for  the 
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first,  third  and  fourth  terms  must  be  zero.  This  leads  to  the  following  matrix  equation 

fill  Cl  0 

-3-113  C2  2 

(2.69) 

9  119  C3  0 

-27  -1  1  27  C4  0 

Solving  (2.69),  the  unknown  coefficients  Ci,C2,C3,  and  C4  are  determined  to  be 
1/24, —9/8,  9/8,  and  —1/24,  respectively.  The  resultant  stencil  of  fourth-order  accuracy 
can  thus  be  written  as 


where  the  (— 3A^^)/(®)(^)  is  the  leading  truncation  error  which  mainly  determines  the 
order  of  accuracy  of  the  stencil.  For  example,  the  fourth-order  stencil  in  the  x  direction  is 
illustrated  in  Fig.  2.7.  Apparently,  in  order  to  evaluate  the  derivative  at  the  center  point, 
one  needs  to  use  the  function  values  of  the  four  neighboring  values.  The  number  above  each 
point  is  the  weighting  coefficient  for  the  summation.  The  method  summarized  in  (2.63)- 
(2.70)  is  referred  to  as  the  method  of  the  undetermined  eoeffieients  [87].  In  general,  this 
method  can  be  used  to  evaluate  the  derivative  not  only  at  the  central  point,  but  also  at  an 


arbitrary  position  in  the  grid. 
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Fig.  2.7.  Fourth-order  central  difference  stencil  along  the  x  direction. 


From  observation,  for  central  difference,  the  only  difference  of  the  coefficients  at  the 
symmetrical  points  is  a  minus  sign.  Therefore,  a  central  difference  stencil  of  A^-th  order  {N 

is  an  even  integer)  can  be  expressed  as  [15],  [86] 

N/2 

f  ^  E  c'n{/[e  +  (2n  -  i)Ae]  -  /[e  -  (2n  -  i)Ae]}  (2.71) 

^  n=l 

Applying  the  method  of  the  undetermined  coefficients  leads  to  the  following  N/2  x  N/2 


linear  system 


1  3 

N-l 

Cl 

1 

1  33  •• 

(A -1)3 

C2 

= 

0 

1  3^-1  •  • 

•  (A -1)^-3 

CnI2 

0 

The  solution  of  (2.72)  can  be  written  in  close  form  as  [86] 


(2.72) 


Cn  = 


(_l)»^-i(iV-  1)!!2 


2^-2(iV/2  +  n-  l)\{N/2  -  n)\{/2n  -  1)2 


(2.73) 


where  x!!  represents  x{x  —  2) (x  —  4)  ... .  The  leading  truncation  error  can  be  written  as 

(2.74) 


Err^  =  f;  C„(2n  - 


2^(N+l)l  ^ 

n=l 
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The  coefficients  and  the  truncation  errors  for  the  central  difference  stencils  up  to  eighth- 
order  accurate  are  summarized  in  Table  2.1  [67],  [86]. 

Table  2.1.  Coefficients  and  Truncation  Errors  of  the  Central  Difference  Stencils 

N  Cl  C-2  Cs  Cj _ Err^ 

2  1  0  0  0 

4  -  ——  0  0  _3AC  fi5)(C) 

^  8  24  ^  ^  640  1  1“’! 

6  75  ^  0  5^f(7)f£) 

^  64  384  640  ^  7168  J 

Q  1225  245  49  _ 5_  35AC 

°  1024  3072  5120  7168  2949124  V?/ 

Notice  that  the  truncation  errors  all  have  the  similar  structure  which  can  be  expressed 
as  [87] 

Err^(A0  ~  (2.75) 

where  C  is  a  constant  determined  by  the  {N  +  l)-th  derivative.  Taking  the  logarithm  of 
both  sides  of  (2.75)  yields 

log  [Err'^ ^  log  C  +  N log  A^  (2.76) 

That  is,  log[Err'^(A^)]  is  a  linear  function  of  logA,^  with  slope  N.  For  example,  let  us 
approximate  the  derivative  of  f{^)  =  cos  (^)  at  ^  =  1.  The  log-log  plots  of  the  truncation 
errors  in  Table  2.1  versus  1/(A^)  are  shown  in  Fig.  2.8.  It  is  evident  that  the  slopes  of  the 
lines  agree  with  the  corresponding  order  of  accuracy.  For  this  case,  increasing  the  order  of 
accuracy  greatly  reduces  the  errors. 
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1/(A^) 

Fig.  2.8.  Truncation  error  of  the  central  difference  stencils. 

To  apply  the  higher-order  scheme,  one  only  need  to  replace  the  second-order  space- 
domain  stencils  in  Yee’s  update  equations  (2.53)  -  (2.58)  by  the  higher-order  stencils  (2.71). 
The  leap-frog  scheme  and  Yee’s  grids  remain  the  same. 

2.4.  Stability  Condition 

In  order  to  guarantee  the  convergence  of  the  numerical  solution,  the  time  step  At 
must  be  bounded.  This  bound  is  usually  referred  to  as  the  Courant  Stability  Condition  [1]. 
Using  the  complex  analysis  of  the  dispersion  relation  [1]  or  the  classical  Von  Neumann 
method  [88],  [14],  the  stability  conditions  for  the  one-,  two-,  and  three-dimensional  Yee 


algorithms  can  be  expressed  as 
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One- dimension: 


A 

At<  — 
c 


(2.77) 


Tw  0- dimensions: 


Three- dimensions: 


Ax2  + 


(2.78) 


(2.79) 


where  At  is  the  time  step,  c  is  the  speed  of  light,  and  Ax,Ay,Az  are  space  steps  with 
respect  to  x,  y,  z  axes,  respectively.  For  the  cubic  mesh  where  (Ax  =  Ay  =  Az  =  As),  the 
stability  conditions  of  (2.77)-(2.79)  reduce  to 


One- dimension: 


As  “ 


(2.80) 


Two- dimensions: 


As  ^/2 


(2.81) 


Three- dimensions: 


S  =  <  Y 

As  ^3 


(2.82) 


where  the  constant  S  is  the  Courant  number  which  describes  the  relation  between  c.  At, 
and  As. 


The  stability  conditions  for  the  higher-order  schemes  have  to  be  derived  individually. 
Here,  the  stability  condition  for  the  important  FDTD(2,4)  scheme  is  summarized  [1] 

One- dimension: 


^  6  Ax 

At  <  -  — 
7  c 


(2.83) 
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Two-dimensions: 


Three- dimensions: 


< - 

/I  +  1 

3 

1 

1 

Aa:^ 


+ 


Ay2 


I 

^  A22 


For  the  cubic  mesh,  (2.83)  -  (2.85)  reduce  to 


One-dimension: 


S 


cAt  6 

-  <  - 

As  -  7 


Tw  0- dimensions: 


cAt  6  1 
~  As  -  7  y/2 


Three- dimensions: 


cAt  6  1 
~  As  -  7 


(2.84) 


(2.85) 


(2.86) 


(2.87) 


(2.88) 


2.5.  Dispersion  Analysis 

The  standard  FDTD  algorithms  are  the  discrete  conterparts  of  Maxwell’s  equations. 
It  is  inevitable  that,  when  the  continuous  derivatives  and  domains  are  discretized,  some 
numerical  errors  are  generated.  Among  these  errors,  the  numerical  dispersion  is  the  most 
critical.  Dispersion  is  a  nonphysical  numerical  effect  meaning  that  the  numerical  phase 
velocity  of  the  wave  is  a  function  of  the  grid  discretization,  operating  frequency,  and  propa¬ 
gation  angles.  As  indicated  in  [1] :  “A  useful  way  to  view  this  phenomenon  is  that  the  FDTD 
algorithm  effectively  embeds  the  electromagnetic  wave  interaction  structure  of  interest  in 
a  tenuous  ‘numerical  aether’  having  a  permittivity  very  close  to  vacuum,  but  not  quite. 


This  ‘aether’  causes  propagating  waves  to  accumulate  delay  or  phase  errors  that  can  lead 
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to  non-physical  results,  such  as  broadening  and  ringing  of  single- pulse  waveforms,  imprecise 
cancelation  of  multiple  scattered  waves,  spurious  anisotropy,  and  pseudo-refraction.”  Thus 
the  accumulation  of  the  dispersion  error  can  be  significant  to  destroy  the  final  results,  es¬ 
pecially  for  problems  that  contain  electrically  large  structures  or  need  very  long  simulation 
time.  Fourier  analysis  has  been  used  to  investigate  the  the  numerical  dispersion  of  a  variety 
of  FDTD  algorithms  [1],  [66],  [67].  In  this  section,  the  dispersion  characteristics  of  the 
standard  FDTD  schemes  will  be  reviewed. 


2.5.1.  The  FDTD(2,2)  Scheme 

Let  us  consider  the  one-dimensional  Yee  algorithm  in  a  free-space  region.  An  ID 
monochromatic  numerical  wave  can  be  expressed  as 


(2.89) 

Hy  =  (2.90) 


where  k  is  the  numerical  wave  number,  K  is  the  space  index  in  the  z-direction  used  to 
distinguish  itself  from  the  wave  number.  Substituting  (2.89)-(2.90)  into  the  finite  difference 
equations  (2.25)-(2.26)  for  the  free-space  case  and  combining  the  resultant  equations  yields 


(2.91) 


Such  a  nonlinear  equation  that  describes  the  relationship  between  the  numerical  wave  num¬ 
ber  k  (or  numerical  phase  speed)  and  the  angular  frequency  to  is  referred  to  as  the  dispersion 
relation.  Apparently,  the  numerical  wave  number  is  no  longer  a  linear  function  of  the  an¬ 
gular  frequency  as  is  the  case  in  the  continuous  world. 
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Mesh  resolution  (N) 


Fig.  2.9.  One-dimensional  dispersion  of  the  Yee  algorithm. 


Using  S  and  N  to  denote  the  Courant  number  cAt/Az  and  the  mesh  resolution  \/Az 
(where  A  is  the  exact  wavelength),  (2.91)  can  be  rewritten  in  a  more  convenient  form  [1],  [66] 


(2.92) 


[w 


=  sm  —  - 
I  N  k 


where  k/k  is  the  ratio  of  the  numerical  wave  number  to  the  exact  wave  number,  which 
describes  the  amount  of  discrepancy  between  them.  Notice  that  k/k  =  1  represents  the 
ideal  no-dispersion  case.  For  one-dimension,  the  dispersion  relation  (2.92)  reduces  to 


k  N  r  1  / 7r5 

v  =  —  sm  —  sm 

k  TT  S  \  N 


(2.93) 


For  different  S  values,  k/k  versus  N  is  plotted  in  Fig.  2.9.  As  shown,  k/k  is  always 
greater  than  unity.  In  other  words,  the  numerical  wave  simulated  by  Yee’s  algorithm  always 
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propagates  slower  than  the  speed  of  light  in  the  continuous  space.  Also  notice  that,  as  N 
increases,  k/k  approaches  unity.  That  is,  the  discrepancy  between  the  numerical  and  exact 
wave  number  becomes  smaller  when  a  finer  mesh  resolution  is  used.  When  S  is  unity 
(Courant  stability  limit),  k  is  always  equal  to  k  for  all  Ns.  This  is  an  ideal  situation 
where  there  is  no  dispersion  error.  The  corresponding  time-step  At  =  Azjc  is  sometimes 
referred  to  as  the  “magic  time-step”  [1].  Unfortunately,  even  if  the  ’’magic  time-step”  is 
used,  the  dispersion  only  vanishes  along  certain  directions  due  to  the  anisotropy  for  multi¬ 
dimensional  cases.  The  anisotropy  of  the  dispersion  can  be  demonstrated  in  2D  and  3D 
dispersion  analyses. 

Following  a  similar  procedure  as  for  the  ID  case,  the  two-dimensional  dispersion  rela¬ 
tion  can  be  obtained  by  substituting  a  monochromatic  traveling-wave  trial  solution  into  the 
2D  Yee’s  algorithms  of  (2.37)-(2.39)  or  (2.40)-(2.42)  and  combining  the  resultant  equations. 


Both  the  TE^  and  TM^  cases  lead  to  the  same  dispersion  relation  of 


1  .  2  /  kxAx  \  1  .  2  /  kyAy 

- — ^  sm  -  4-  - — ^  sm  — - 

I  2  i  Ay2  I  2 


(2.94) 


where 


kx  =  k  cos  (j),  ky  =  k  sin  cj) 


Notice  that,  in  two  dimensions,  the  dispersion  relation  is  also  a  function  of  the  propagation 
angle  (j).  For  the  square  mesh  (Ax  =  Ay  =  As),  the  two-dimensional  dispersion  relation 
can  be  rewritten  as 


1  9  /  7r5  \  9  /  vr  cos  d)k\  9  /  vr  sin  6  k 

=.m 


(2.95) 


The  two-dimensional  dispersion  relation  in  (2.95)  is  a  transcendental  function  which  has  no 


close  form  solution  for  k/k.  However,  since  the  derivative  of  the  transcendental  function  is 
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readily  obtained,  (2.95)  can  be  numerically  solved  using  the  Newton-Raphson’s  method  [89]. 
The  procedure  is  summarized  below. 

Letting  x  =  fc//c,  a  =  ncoscp/N,  b  =  vrsin  and  d  =  sin^(7r5/A^)/5^,  (2.95)  can  be 
rewritten  as 

f{x)  =  sin^(ax)  +  sin^(6x)  —  =  0  (2.96) 

The  root  x  of  (2.96)  can  be  solved  by  performing  the  following  iteration 

fi^i)  (r,  orTN 

Xi+i  =  Xi-  (2.97) 

f  [Xi) 

where  f'{x)  is  the  derivative  of  f{x)  which  takes  the  form  of 

f\x)  =  asin(2ax)  +  6sin(26x)  (2.98) 

Performing  the  iteration  of  (2.97),  k/k  can  be  readily  solved. 

Setting  S  =  0.5,  k/ks  versus  the  propagation  angle  0  for  various  mesh  resolutions  N 
is  plotted  in  Fig.  2.10.  Since  the  dispersion  is  90°  periodic,  only  the  curves  within  the  first 
quadrant  are  plotted.  It  is  evident  that  all  the  dispersion  curves  are  anisotropic  (functions 
of  propagation  angles)  which  possess  their  maxima  along  the  principal  axes  0°  and  90°  and 
the  minima  along  the  diagonal  axes  45°.  Similar  to  the  one-dimensional  case,  the  dispersion 
decreases  as  the  mesh  resolution  N  increases. 

When  the  mesh  resolution  is  fixed  at  N  =  10,  k/k  versus  the  propagation  angle  (j) 
for  various  S  are  plotted  in  Fig.  2.11.  As  S  approaches  the  2D  stability  limit  (l/\/2), 
the  dispersion  curve  shifts  down  in  parallel  towards  unity.  In  other  words,  all  the  curves 
maintain  the  same  shape  for  different  S.  When  3  =  1/ \/2,  the  dispersion  eventually  vanishes 


along  the  diagonal  axis  (j)  =  45°. 


Propagation  angles  (p) 


Fig.  2.10.  Two-dimensional  dispersion  of  the  Yee  algorithm  (5  =  0.5) 


Fig.  2.11.  Two-dimensional  dispersion  of  the  Yee  algorithm  (N  =  10) 
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The  three-dimensional  dispersion  relation  of  Yee’s  algorithm  can  be  readily  derived 
and  written  as 


ky/\y 


kz^z 


(2.99) 


where 


kx  =  kcos  (j)sm9,  ky  =  k  sin  4>  sin  9,  kz  =  k  cos  9, 


For  the  cubic  mesh  (Ax  =  Ay  =  Az  =  As),  (2.99)  can  be  rewritten  in  the  form  of 


2  /  vr  cos  (j)  sin  9  k\  2  /  sin  (j)  sin  9  k\  2 
- W - 1  - N - 


TT  COS  9  k 
N  k 


(2.100) 


For  a  fixed  Courant  number  S  =  0.4,  k/k  versus  propagation  angles  {9,(p)  for  N  =  5 
and  =  10  is  plotted  in  Figs.  2.12  and  2.13,  respectively.  As  N  increases,  the  absolute 
value  and  the  anisotropy  of  the  dispersion  are  significantly  reduced. 

For  a  fixed  mesh  resolution  N  =  10,  k/ks  versus  propagation  angles  (0,  (p)  for  different 
Courant  numbers  (5  =  0.1  and  S  =  0.577)  is  plotted  in  Figs.  2.14  and  2.15,  respectively. 
It  is  evident  that  the  shapes  of  the  dispersion  for  different  S  are  identical.  As  S  approaches 
the  Courant  limit  (l/\/3),  the  k/k  distribution  shifts  down  towards  unity. 
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Fig.  2.12.  Three-dimensional  dispersion  of  the  Yee  algorithm  (5  =  0.4,  =  5). 


Fig.  2.13.  Three-dimensional  dispersion  of  the  Yee  algorithm  (S'  =  0.4,  Ai  =  10). 
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0  <1’ 


Fig.  2.14.  Three-dimensional  dispersion  of  the  Yee  algorithm  (S'  =  0.1,  =  10). 


0  <1’ 

Fig.  2.15.  Three-dimensional  dispersion  of  the  Yee  algorithm  (S  =  0.577,  Y  =  10). 
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2.5.2.  The  FDTD(2,4)  Scheme 


For  conciseness,  only  the  three-dimensional  dispersion  of  FDTD(2,4)  is  discussed  in 
this  section.  The  3D  dispersion  relation  for  the  FDTD(2,4)  scheme  can  be  written  as 


Sin 


("T) 


1 


Ax^ 

1 

1 


9  . 

-  sm 


-  sm 


kxAx 


kyAy 


1  .  ( 3kxAx 

- sm  - 

24  I  2 
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where 


kx  =  k  cos  (f>  sin  0 ,  ky  =  k  sin  (j)  sin  0 ,  kz  =  k  cos  9, 


For  the  cubic  mesh  (Ax  =  Ay  =  Az  =  As),  (2.101)  can  be  rewritten  in  the  form  of 
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(2.102) 


Notice  that  the  two-dimensional  dispersion  relation  can  be  readily  obtained  by  forcing 
9  =  90°  in  (2.101)  or  (2.102).  The  one-dimensional  dispersion  relation  can  be  obtained  by 
forcing  {9  =  90°,  cj)  =  0°)  in  (2.101)  or  (2.102).  The  one-dimensional  dispersions  versus  mesh 
resolution  N  for  different  S  are  plotted  in  Fig.  2.16.  As  shown,  all  the  curves  approach 
unity  when  N  increases.  However,  the  dispersion  characteristic  of  the  FDTD(2,4)  scheme 
is  very  different  from  that  of  the  FDTD  (2,2)  scheme.  There  is  an  optimum  S  value  {Sopt) 
for  each  N  where  the  dispersion  is  minimized  to  zero  [66].  When  S  deviates  from  Sopt,  the 
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Fig.  2.16.  One-dimensional  dispersion  of  the  FDTD(2,4)  algorithm. 
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Fig.  2.17.  Optimal  S  for  the  ID  FDTD(2,4)  algorithm. 
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discrepancy  between  k  and  k  increases.  Also  notice  that  if  S'  <  Sopt,  k  >  k]  \i  S  >  Sopt, 
k  <  k.  The  optimum  S  value  can  be  found  by  forcing  k  =  k  m  the  dispersion  relation  and 
solving  the  resultant  nonlinear  equation.  For  the  ID  FDTD(2,4),  Sopt  versus  N  is  plotted  in 
Fig.  2.17.  Notice  that  the  expected  fourth-order  convergence  is  not  observed  in  Fig.  2.16. 
The  convergence  rate  of  the  FDTD(2,4)  is  only  second-order.  This  can  be  attributed  to  the 
fact  that  the  error  generated  by  the  time-domain  discretization  is  still  second  order,  which 
deteriorates  the  accuracy  of  the  entire  scheme  unless  At  is  sufficiently  small. 

Next,  the  2D  dispersion  of  the  FDTD(2,4)  scheme  will  be  discussed.  For  S  =  0.5,  the 
2D  dispersion  versus  propagation  angles  for  different  N  is  plotted  in  Fig.  2.18.  Different 
from  Yee’s  algorithm,  the  maximum  discrepancy  occurs  along  45°,  and  the  minima  are  along 
the  principal  axes.  As  N  increases,  the  absolute  value  of  the  dispersion  and  the  anisotropy 
decreases.  For  a  fixed  mesh  resolution  N  =  10,  the  dispersion  versus  propagation  angles  for 
different  S  is  plotted  in  Fig.  2.19.  Similar  as  the  FDTD(2,2)  case,  the  shapes  of  all  curves 
remain  identical.  When  S  deviates  the  optimum  value,  the  dispersion  curves  shifts  away 
from  unity. 

To  complete  the  discussion,  the  3D  dispersion  plots  of  the  standard  FDTD(2,4)  scheme 
are  also  presented.  Firstly,  the  Courant  number  is  fixed  to  S  =  OA,  k/k  versus  propagation 
angles  for  N  =  5  and  A  =  10  is  plotted  in  Figs.  2.20  and  2.21,  respectively.  It  is  evident 
that  using  a  finer  mesh  resolution  greatly  reduces  the  anisotropy  of  the  dispersion.  Secondly, 
when  the  mesh  resolution  is  fixed  to  N  =  10,  k/k  versus  propagation  angles  for  S  =  0.49 
and  S  =  0.15  is  plotted  in  Figs.  2.22  and  2.23,  respectively.  As  shown,  the  shapes  of 


dispersion  are  identical  but  they  are  shifted. 
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Fig.  2.18.  Two-dimensional  dispersion  of  Fang’s  algorithm  (5  =  0.5) 


Fig.  2.19.  Two-dimensional  dispersion  of  Fang’s  algorithm  (N  =  10) 
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0  <1’ 


Fig.  2.20.  Three-dimensional  dispersion  of  the  FDTD(2,4)  algorithm  (N  =  5). 


0  <1’ 

Fig.  2.21.  Two-dimensional  dispersion  of  the  FDTD(2,4)  algorithm  (N  =  10). 
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Fig.  2.22.  Three-dimensional  dispersion  of  the  FDTD(2,4)  algorithm  (S  =  0.49). 


Fig.  2.23.  Two-dimensional  dispersion  of  the  FDTD(2,4)  algorithm  (5  =  0.15). 
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2.5.3.  The  FDTD(2,N)  Schemes 


The  dispersion  relation  for  any  FDTD(2,N)  scheme  can  be  expressed  as  [86] 
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(2.103) 


The  left-hand-side  of  (2.103)  depicts  the  time-domain  discretization  error  while  the  right- 
hand-side  describes  the  error  generated  by  the  space-domain  discretization.  For  the  stan¬ 
dard  FDTD(2,N)  algorithms,  three  important  arguments  can  be  made: 


•  Both  the  space  and  time  discretization  errors  contribute  to  the  total  dispersion  of  the 
scheme.  In  order  to  achieve  the  highest  possible  accuracy,  it  is  important  to  balance 
the  orders  of  accuracy  of  the  time-  and  space-discretizations. 


•  The  time-domain  discretization  errors  are  identical  along  any  propagation  direction 
(isotropic).  The  effect  of  varying  the  Courant  number  is  to  shift  the  dispersion  curve. 
Depending  on  the  schemes,  there  is  an  optimum  S  value  that  minimize  the  dispersion 
for  each  mesh  resolution. 


•  The  errors  generated  by  the  space  discretization  are  anisotropic.  The  use  of  a  finer 
mesh  resolution  decreases  the  anisotropy. 

In  order  to  achieve  a  high  accuracy,  ideally  the  time  domain  scheme  needs  to  be  at  least 
N-th  order.  However,  for  simplicity,  a  popular  treatment  is  to  use  an  FDTD(2,N)  scheme 
with  a  small  time-step  [33].  Moreover,  the  space  domain  scheme  is  more  critical  since  it 
controls  both  the  absolute  value  and  the  anisotropy  of  the  dispersion. 

In  summary,  the  derivation  of  the  standard  FDTD  schemes  are  based  on  the  Tay¬ 


lor  series  expansion.  The  coefficients  of  the  stencils  are  selected  to  minimize  the  leading 
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truncation  error.  These  coefficients  might  not  be  the  optimal  in  terms  of  the  dispersion 
error.  Alternative  methods,  with  the  ultimate  goal  of  minimizing  the  dispersion  error  are 
still  highly  desirable.  Such  methods  include  the  nonstandard  FDTD  schemes  which  are  con¬ 
structed  directly  based  on  the  dispersion  analysis.  In  the  following  chapter,  the  fundamental 
theory  of  the  nonstandard  FDTD  schemes  will  be  introduced. 


CHAPTER  3 


THE  NONSTANDARD  EINITE-DIEEERENCE  METHODS 

In  Chapter  2,  we  introduced  the  fundamental  theory  of  the  Standard  Einite-Difference 
Methods  (SED),  from  the  simplest  Yee’s  algorithm  to  the  more  complex  Higher-Order 
EDTD  (HOEDTD).  Using  the  Taylor  series  expansion,  the  coefficients  for  the  SED  stencils 
are  derived  by  minimizing  the  leading  truncation  error.  Usually,  the  advantage  of  the 
HOEDTD  becomes  obvious  only  if  the  mesh  resolution  is  refined.  However,  for  electrically 
large  problems,  a  very  coarse  mesh  resolution  usually  has  to  be  used  due  to  the  tight  budget 
of  the  available  computer  memories.  In  addition,  for  the  unbalanced  schemes,  such  as  the 
EDTD(2,4)  scheme,  the  gain  of  accuracy  in  the  space  domain  could  be  deteriorated  by  the 
large  time-domain  errors  unless  a  sufficiently  small  time  step  is  used.  Consequently,  it  is 
questionable  whether  or  not  the  coefficients  for  the  SED  stencils  are  optimal  in  terms  of 
dispersion.  Instead  of  the  pursuit  of  the  higher-order,  alternative  methods  with  the  ultimate 
goal  of  low-dispersion  are  more  practical.  The  Nonstandard  Einite-Difference  Methods 
(NSED)  are  such  schemes  that  are  constructed  directly  based  upon  the  dispersion  analysis. 
The  concept  of  the  NSED  was  first  proposed  by  Mickens  [37].  Then,  Cole  applied  the  NSED 
to  solve  the  Maxwell’s  equations  [40],  [41].  In  this  report,  the  NSED  methods  refer,  in  its 
broadest  sense,  to  those  schemes  which  are  modified  based  upon  the  SED  schemes  to  further 
reduce  the  dispersion.  To  distinguish  with  the  SED  methods,  an  NSED  scheme  using  M- 
and  N-point  stencils  in  the  time  and  space  domains  will  be  designated  as  the  (NS)MN 
scheme.  Notice,  (M,  N)  no  longer  represents  the  order  of  accuracy. 

In  this  chapter,  the  fundamentals  of  the  NSED  method  will  be  first  reviewed  for  the 


ID  case.  Then  the  performance  of  the  basic  ID  NS22  scheme  and  its  direct  extension,  the 
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NS24  scheme,  will  be  discussed.  For  multi-dimensional  cases,  an  improved  NSFD  method 
(INS)  will  be  presented.  Furthermore,  the  concept  of  the  NSFD  will  be  generalized  to 
construct  a  more  General  NSFD  method  (GNS).  Both  the  dispersion  analysis  and  the  nu¬ 
merical  simulations  will  be  demonstrated.  Afterwards,  Gauss’s  Law  for  the  NSFD  methods 
will  be  examined.  Finally,  to  reduce  the  errors  generated  at  the  material  discontinuities, 
two  material  interface  conditions  will  be  developed  for  the  NSFD  methods  with  extended 
stencils. 


3.1.  The  One-Dimensional  NS22  Scheme 

For  clarity,  the  discussion  of  the  NSFD  methods  starts  with  the  one-dimensional  NS22 
scheme,  which  uses  a  two-point  stencil  in  both  the  time  and  space  domains. 


3.1.1.  Fundamental  of  the  NS22  scheme 

The  standard  way  to  approximate  a  partial  derivative  is  to  use  the  second-order  central- 
difference  operator  introduced  in  Ghapter  2.  For  example,  the  space  derivative  with  respect 
to  X  can  be  written  as 
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where  represents  the  numerator  of  the  standard  second-order  central-difference.  If 
one  substitutes  the  monochromatic  wave-equation  solution  /o  =  into  it  turns 
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2  2  .  /  /cqAx 

Clearly,  the  result  on  the  right  hand  side  is  the  exact  solution  multiplied  by  a  constant  sk'^Q 
Dividing  both  sides  of  (3.2)  by  yields 
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That  is,  the  new  operator  =  [(d^”'^)/(sA:^o)]  gives  the  exact  solution  at  one  single 

frequency  (ko)  along  the  xaxis.  Based  on  this  observation,  the  space- domain  NS22  operators 
are  defined  as 
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A  similar  procedure  is  followed  to  derive  the  time-domain  NS22  operator  which  can  be 
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Notice  that  sA:|q  depends  on  ko  and  =  x,  y,  z);  suq  is  determined  by  ujo  and  At.  Thus, 
the  correction  parameters  s/c|q  and  suJq  can  be  considered,  respectively,  as  the  frequency- 
optimized  space  and  time  increments.  To  apply  the  NS22  algorithm,  one  only  needs  to 
replace  the  space  and  time  steps  (A^,  At)  in  Yee’s  algorithm  by  their  frequency-optimized 
counterparts  s/cq  and  swg,  respectively. 
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3.1.2.  Dispersion  Analysis 

In  this  section,  we  will  explore  the  basic  properties  of  the  ID  NS22  scheme  by  analyzing 
its  dispersion  characteristics.  Following  a  similar  procedure  introduced  in  Section  2.5,  the 
ID  dispersion  relation  for  the  NS22  scheme  can  be  derived  and  written  as 


where  u)q  is  the  design  angular  frequency  and  /cq  is  the  corresponding  exact  wave  number 
at  the  design  frequency.  Equation  (3.7)  can  be  rewritten  as 

(^)  it)  ^  (f ) 

where  S  =  cAz/At  is  the  Courant  number;  N  and  Nq  represent  the  mesh  resolution  at  an 
arbitrary  and  at  a  design  frequency,  respectively.  Notice  that  when  N  =  Nq,  the  two  terms 
containing  S  in  (3.7)  cancel  out,  which  imply  zero  dispersion  {k/k  =  1)  at  N  =  Nq. 

To  better  demonstrate  the  characteristic  of  the  dispersion,  the  local  dispersion  error  is 
defined  as 

Eioc  =  (^1  - 

Fig.  3.1  illustrates  the  absolute  values  of  the  local  dispersion  errors,  for  a  fixed  design  mesh 
resolution  Nq  =  10  and  various  Courant  numbers  S.  The  local  dispersion  errors  of  Yee’s 
algorithm  are  also  plotted  as  references.  The  figure  clearly  shows  that,  for  the  NS22  scheme, 
the  dispersion  error  is  extremely  low  at  N  =  Nq.  However,  compared  to  the  errors  generated 
by  Yee’s  algorithm,  this  excellent  low-dispersion  performance  only  presents  in  the  vicinity 
of  Nq.  If  N  is  increased  from  Nq  to  infinity,  the  dispersion  error  of  the  NS22  scheme  could 


be  even  worse  than  that  of  Yee’s  algorithm.  Similar  to  the  Yee  algorithm,  using  a  Courant 
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Fig.  3.1.  One-dimensional  local  dispersion  error  of  the  NS22  scheme  {Nq  =  10). 


number  which  is  as  close  as  possible  to  the  Courant  limit  helps  reduce  the  overall  dispersion 
error  of  the  NS22  scheme.  Therefore,  the  NS22  scheme  is  a  narrow-band  algorithm  in  terms 
of  the  dispersion,  which  is  suitable  for  single-frequency  simulations.  The  local  dispersion 
errors  versus  N  for  a  fixed  Courant  number  S  =  0.5  and  Nq  =  10,  20  are  plotted  in  Fig.  3.2. 
It  is  evident  that  the  dispersion  null  of  the  NS22  scheme  can  be  readily  shifted  by  choosing 
different  Nq.  As  Nq  approaches  infinity,  the  NS22  curve  asymptotically  approaches  the  curve 
of  Yee’s  algorithm.  Notice  that,  for  N  <  Nq  (higher  frequency  region  for  a  fixed  space  step), 
the  dispersion  error  of  the  ID  NS22  scheme  is  always  smaller  than  that  of  the  Yee  algorithm. 
Bearing  in  mind  that  the  higher-frequency  signals  usually  generate  greater  dispersion  error. 
This  fact  indicates  that  the  NS22  scheme  still  can  be  used  for  the  broad-band  simulations 


if  one  selects  an  appropriate  Nq. 
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Fig.  3.2.  One-dimensional  local  dispersion  error  of  the  NS22  scheme  (5  =  0.5). 

3.1.3.  Numerical  Simulations 

In  this  section,  the  characteristics  of  the  ID  NS22  scheme  will  be  demonstrated  through 
numerical  simulations.  As  a  first  example,  a  continuous  sinusoid  wave  with  /o  =  2GHz  is 
generated  at  the  center  of  an  ID  free-space  region.  The  ID  space  is  large  enough  so  that 
the  traveling  waves  will  never  reach  the  domain  boundaries  within  the  duration  of  interest. 
The  space  step  and  the  Courant  number  are  set  to  Az  =  0.01  m  (A/15  at  /o)  and  S  =  0.5, 
respectively.  A  snapshot  is  taken  for  the  magnitude  of  the  E  field  at  t  =  lOOOAt.  The 
snapshots  for  both  the  Yee  and  the  NS22  schemes  are  zoomed  near  2OOA2;  and  plotted  in 
Fig.  3.3.  It  is  evident  that  the  solution  predicted  by  the  NS22  has  an  excellent  agreement 
with  the  analytical  solution.  However,  the  results  provided  by  Yee’s  algorithm  exhibits  a 
phase  lag  of  about  24°.  For  a  longer  simulation  time,  this  phase  discrepancy  is  even  greater 
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Fig.  3.3.  Propagation  of  the  ID  sinusoid  waves  {Az  =  0.01m, 5  =  0.5). 

due  to  the  accumulation  of  the  dispersion.  However,  if  5  =  1,  the  ID  Yee’s  algorithm  and 
NS22  scheme  both  produce  the  exact  solution. 

A  second  simulation  is  a  half-wavelength  dielectric  slab  (e^  =  16)  in  a  free-space  region 
[37].  The  space  step  and  the  Courant  number  are  set  to  Az  =  0.0025m  and  S  =  0.99  in 
terms  of  the  speed  of  light  in  free-space  region,  respectively.  The  incident  wave  is  a  sequence 
of  a  sinusoid  wave  with  /o  =  3GHz,  as  shown  in  the  first  subplot  of  Fig.  3.4.  Ideally,  there  is 
no  reflection  if  the  sequence  is  infinitely  long.  An  average  permittivity  is  applied  along  the 
material  interface.  Notice  that  there  is  almost  no  dispersion  in  the  free-space  region  since  S 
is  very  close  to  the  “magic  time  step”  S'  =  1.  The  errors  are  dominated  by  the  phase  error 
inside  the  slab  and  the  nonphysical  errors  generated  by  the  material  discontinuities.  The 
reflected  and  the  transmitted  E  fields  at  t  =  llOOAt,  predicted  by  the  Yee  and  the  NS22 
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schemes,  are  plotted  in  the  second  and  third  subplots  of  Fig.  3.4.  As  shown,  large  reflections 
are  observed  at  both  ends  of  the  reflected  sequence  due  to  the  signal  discontinuities  there. 
As  shown,  the  results  predicted  by  Yee’s  algorithm  exhibit  large  reflection  along  the  center 
part  of  the  reflected  wave.  This  can  be  attributed  by  the  dispersion  of  Yee’s  algorithm  which 
makes  the  thickness  of  the  slab  no  longer  half  wavelength.  However,  the  NS22  scheme  has 
the  freedom  to  choose  the  sk  value  inside  the  material  which  automatically  adjust  the 
phase  speed.  It  is  evident,  the  results  predicted  by  the  NS22  scheme  present  much  smaller 
reflections.  Based  on  the  above  observations,  NSFD  algorithm  is  suitable  for  problems 
which  require  accurate  phase  cancelation  at  a  single  operating  frequency. 

Both  of  the  previous  examples  are  single-frequency  simulations.  However,  whether  the 
NSFD  algorithm  is  capable  of  simulating  broad-band  signals  remains  a  question.  To  explore 
the  broad-band  performance  of  the  NSFD  algorithm,  a  Gaussian  pulse  with  (a  =  30,  /3  =  10) 
is  generated  in  the  ID  free-space  region.  The  space  step  and  the  Courant  number  are  set 
to  Az  =  0.01m  and  S  =  0.5,  respectively.  The  frequency  spectrum  of  the  Gaussian  pulse 
is  plotted  in  Fig.  3.5.  As  shown,  the  Gaussian  pulse  contains  significant  low-frequency 
components.  The  -3dB  frequency  of  the  Gaussian  pulse  is  located  at  about  f'j.dB  =  1-6 
GHz.  For  the  NS22  scheme,  the  design  frequency  is  set  to  fi,dB-  The  snapshots  of  the  E 
field  at  t  =  500At  for  the  Yee  and  the  NS22  schemes  are  plotted  in  Fig.  3.6.  As  shown,  the 
pulses  predicted  by  the  Yee  and  NS22  schemes  suffer  from  similar  level  of  phase  errors.  This 
example  suggests  that,  by  properly  choosing  a  design  mesh  resolution,  the  NSFD  algorithm 


may  generate  dispersion  errors  no  worse  than  Yee’s  algorithm. 


Magnitude  of  the  E  Field 


Frequency  (GHz) 


Fig.  3.5.  Frequency  spectrum  of  a  Gaussian  pulse  (a  =  30, /3  =  10). 
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Fig.  3.6.  Propagation  of  the  ID  Gaussian  pulse  (Az  =  0.01m, 5  =  0.5) 
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3.2.  The  One-Dimensional  NS24  Scheme 


The  basic  concept  of  the  NS22  scheme  can  be  readily  extended  to  Fang’s  FDTD  (2,4) 
algorithm  [15]  to  construct  a  NS24  scheme  [90],  [91].  Following  a  similar  procedure,  the 
space-domain  NS24  operator  can  be  defined  as 
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Similar  to  the  NS22  case,  s/c|q  can  be  considered  as  the  frequency-optimized  space  step. 
Since  the  same  second-order  stencil  is  used  in  the  time  domain,  the  time-domain  NS24 
operator  remains  the  same  as  in  (3.5).  To  apply  the  NS24  algorithm,  one  needs  to  replace 
the  space  and  time  steps  (A^,  At)  in  the  regular  FDTD  (2,4)  algorithm  by  their  frequency- 
optimized  counterparts  sk^  and  suJq,  respectively. 

The  dispersion  relation  of  the  ID  NS24  scheme  can  be  readily  derived  and  expressed 


as 


1  f  SkoAz 

- sm  - 

24  V  2 


3kAz 

2 


2 


(3.12) 


Comparing  to  the  dispersion  relation  of  the  ID  NS22  scheme  in  (3.7),  the  only  difference 
is  the  first  term  on  the  left  hand  side.  Similarly,  the  dispersion  relation  of  (3.12)  can  be 
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rewritten  as 


(3.13) 


It  can  be  expected  that  the  dispersion  error  curve  of  the  NS24  scheme  will  present  a  null 
which  is  located  at  =  Nq.  It  has  been  indicated  in  Section  2.5  that,  for  the  FDTD(2,4) 
scheme,  there  could  be  a  second  null  which  is  controlled  by  the  Courant  number  S.  To 
control  the  position  of  the  two  nulls  of  the  NS24  scheme,  we  can  specify  two  design  mesh 
resolutions  Nq  and  Ni.  The  desired  S  value  can  be  found  as  follows.  Firstly,  k/k  is  forced 
to  unity  in  (3.13).  Then  substituting  N  =  Ni  into  (3.13)  yields 


(3.14) 


Solving  (3.14)  using  the  Newton-Ralphson  method,  the  desired  S  value  can  be  readily 
obtained.  Notice  that,  due  to  the  symmetry  of  the  equation,  switching  Nq  and  A^i  in  (3.14) 
will  not  affect  the  result.  The  computed  S  values  for  various  {Nq  and  Ni)  combinations  are 
summarized  in  Table  3.1. 


Table  3.1.  Computed  S  values  for  various  Nq  and  Ni 


Nq 

Ni 

S 

Case  1 

17.5 

17.5 

0.1693 

Case  2 

15 

20 

0.1745 

Case  3 

25 

10 

0.2243 

The  local  dispersion  errors  of  each  of  the  Case  in  Table  3.1  are  plotted  in  Figs.  3.7  -  3.9, 
respectively.  It  is  evident  that  the  dispersion  nulls  are  located  precisely  at  the  design  mesh 


Fig.  3.7.  One-dimensional  local  dispersion  error  of  the  NS24  scheme  {Nq 
0.1693). 


17.5,5  = 


Fig.  3.8.  One-dimensional  local  dispersion  error  of  the  NS24  scheme  {Nq  =  15,  5  =  0.1745) 
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Fig.  3.9.  One-dimensional  local  dispersion  error  of  the  NS24  scheme  {Nq  =  25,  S  =  0.2243). 

resolutions  {Nq  and  A^i).  By  choosing  different  combinations  of  A^o  and  it  is  possible 
to  control  the  low-dispersion  band  width  for  the  ID  NS24  scheme. 

3.3.  The  Anisotropy  for  the  Multi-Dimensional  NSFD  schemes 

The  NSFD  algorithms  discussed  so  far  are  limited  to  the  one-dimensional  domain. 
For  the  multi-dimensional  cases,  the  ideal  performance  of  the  ID  NSFD  schemes  could  be 
totally  deteriorated  by  another  nonphysical  artifact:  the  anisotropy.  To  best  demonstrate 
this  issue,  we  start  the  discussion  with  the  dispersion  analysis  for  the  multi-dimensional 
NS22  scheme. 

The  multi-dimensional  NSFD  schemes  are  constructed  by  simply  replacing  the  time 
and  space  increments  (At  and  A^)  by  their  nonstandard  counterparts.  Following  a  similar 
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procedure,  the  dispersion  relations  of  the  multi-dimensional  NSFD  schemes  can  be  readily 
obtained.  The  general  3D  dispersion  relation  of  the  NS22  scheme  can  be  expressed  as 

where 

kx  =  kcoscpsmO,  ky  =  k  sin  cj)  sin  6 ,  kz  =  k  cos  6. 

For  the  cubic  mesh  case  (Ax  =  Ay  =  Az  =  As),  this  relation  reduces  to 

it)  =  (^t) 

where  k/kis  the  ratio  of  the  numerical  wave  number  to  the  exact  wave  number;  Nq  =  Aq/As 
is  the  design  mesh  resolution.  The  2D  dispersion  relations  can  be  readily  obtained  by  forcing 
(^  =  0°  in  (3.16).  Clearly,  the  dispersion  error  for  the  multi-dimensional  NSFD  schemes  is 
also  a  function  of  (j). 

The  2D  local  dispersion  errors  versus  the  (j)  for  the  2D  Yee  and  NS22  schemes  are 
plotted  in  Fig.  3.10.  Only  the  first  quarter  is  plotted  since  all  the  curves  are  90°  periodic. 
Although,  the  NS22  scheme  exhibits  zero  dispersion  along  the  principle  axes  (0°  and  90°), 
it  still  suffers  from  a  similar  level  of  anisotropy  as  does  the  Yee  algorithm. 

3.4.  The  Improved  NSFD  Schemes 

In  order  to  reduce  the  average  dispersion  error,  a  family  of  methods  [1] ,  [53]-  [55] ,  [60] 
was  constructed  based  on  shifting  the  symmetrical  center  of  the  dispersion  error  curve  to 
zero  at  a  design  mesh  resolution.  Using  a  similar  approach,  an  improved  NSFD  method 
(INS)  [91]  is  presented  in  this  section.  By  introducing  extra  degrees  of  freedom  into  the 
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Fig.  3.10.  Two-dimensional  local  dispersion  error  of  the  NS22  scheme  {N  =  Nq  =  10,5  = 
0.7). 

spatial  operators  of  the  regular  NSFD  methods  (the  NS22  and  INS24  schemes),  the  INS22 
and  INS24  methods  are  able  to  mitigate  the  average  dispersion  error. 

3.4.1.  The  INS22  and  INS24  Schemes 

Since  the  derivation  of  the  regular  NSFD  schemes  only  considers  the  plane-wave  prop¬ 
agating  along  the  principle  axes,  sA:|q  and  sA:|q  may  not  be  the  optimal  space  steps  for  waves 
traveling  in  other  directions.  To  achieve  the  optimal  correction  parameters  in  the  global 
sense,  one  can  introduce  extra  degrees  of  freedom  into  (3.4)  and  define  the  space-domain 
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INS22  operators  as 


j^ins22  f  _  /ir+1/2  fx  1/2 


j~\ins22  p  _  4+1/2  -  4-1/2 
r^ins22  f  _  4+1/2  4—1/2 


(3.17) 

(3.18) 


(3.19) 


where  a,  b,  and  c  are  extra  free  parameters.  For  simplicity,  we  only  consider  the  cubic  mesh 
where  Ax  =  Ay  =  Az  =  As.  In  that  case,  the  free  parameters  are  identical  (a  =  6  =  c  =  1). 
Thus  the  space  domain  INS22  operators  in  (3.17)  -  (3.19)  reduce  to 


^ins22f  _  4+1/2  -  4-1/2 


r\inszz  £ 

4  = 


q  ■  skL 


,  i  =  x,y,z 


(3.20) 


Following  a  similar  procedure,  the  space-domain  INS24  operators  can  be  expressed  as 


T^ins2if  _  I  (4+1/2  -  4-1/2)  -  M  (4+3/2  “  4-3/2)  ^  ^ 


(3.21) 


Notice  that  the  regular  NSFD  schemes  are  special  cases  of  the  improved  NSFD  schemes 
when  the  free  parameter  ^  =  1.  The  values  of  q  can  be  determined  based  on  the  dispersion 
analysis  that  follows. 


3.4.2.  Dispersion  Analysis  for  INS22 

The  two-dimensional  dispersion  relation  of  the  INS22  scheme  can  be  found  by  substi¬ 
tuting  a  harmonic  solution  to  the  2D  INS22  update  equations.  Following  a  similar  notation 
as  in  [67],  the  2D  INS22  dispersion  relation  can  be  written  as 


sin^(7r/Ao)  o  f  7rS\  1  o  /  7rcos(/>/c\  n 

— TT-, ^ — r  sm  — —  =  ^  sm  — — —  -  M-  sm 

sin^(7rS'/Ao)  J  \  N  k  j 


TT  sin  (p  k 
N  k 


(3.22) 
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where  k  is  the  numerical  wave  number,  k  is  the  exact  wave  number,  N  represents  the  actual 
cell  numbers  per  wavelength,  Nq  represents  the  cell  numbers  per  wavelength  at  the  design 
frequency,  and  S  =  (cAt)/As  is  the  Courant  number. 

Ideally  at  the  design  mesh  resolution  Nq,  we  desire  the  dispersion  error  to  be  equal  to 
zero  for  all  propagation  angles.  That  is,  in  (3.22),  k  should  be  equal  to  k  for  all  angles  in 
the  range  of  (/>  =  0°  ~  90°  (note  that  the  dispersion  curve  is  90°  periodic)  at  N  =  Nq.  For 
a  finite  number  of  angles  (/),  this  leads  to  an  overdetermined  linear  system  (the  number 
of  the  equations  is  greater  than  the  number  of  the  unknowns).  Assuming  the  angles  are 
uniformly  sampled  with  A(p  =  7r/[2(/  —  1)],  the  resultant  I  equations  can  be  expressed  as 


(3.25) 


(3.26) 
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Fig.  3.11.  Errors  of  q  for  the  2D  INS22  scheme. 


Notice  that  when  N  =  Nq,  the  two  terms  containing  S  in  the  left-hand-side  of  (3.22)  cancel 
out  which  makes  q  independent  of  At. 

A  remaining  question  is  how  many  angles  we  need  to  calculate  ql  Using  the  q  value 
with  large  I  {I  =  90,  001)  as  a  reference,  the  errors  of  q  versus  different  I  are  plotted  in  Fig. 
3.11.  As  shown,  the  error  of  q  decreases  as  I  or  Nq  increases.  For  Nq  =  10,  to  achieve  an 
error  for  q  smaller  than  10“®,  I  must  be  larger  than  400  (/  >  400). 

The  local  dispersion  errors  of  the  three  (2,2)  algorithms  are  plotted  in  Fig.  3.12.  Only 
the  curves  between  0°  and  90°  are  plotted  because  all  the  curves  are  90°  periodic.  As  shown, 
the  average  dispersion  error  of  NS22  is  smaller  than  that  of  Yee’s  algorithm.  Also,  the  2D 
NS22  scheme  has  zero  dispersion  along  0°  and  90°.  For  the  2D  INS22  scheme,  the  dispersion 
vanishes  at  22.5°  and  67.5°  due  to  the  symmetry.  The  effect  of  the  free  parameter  q  in  the 


INS22  scheme  is  to  shift  the  average  error  down  to  zero.  However,  the  maximum  phase 
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Fig.  3.12.  Two-dimensional  local  dispersion  errors  of  the  (2,2)  schemes  [S  =  0.6,  N  =  Nq  = 
10, q{I  =  401)  =  1.00414816]. 


errors  (peak-to-peak)  of  the  above  three  schemes  are  the  same. 

Similar  to  [1],  [53]-  [55],  [60],  forcing  k  =  k  at  N  =  Nq  along  a  zero  dispersion  angle 
{4>o)  in  (3.22)  leads  to 


q  = 


\/ sin^(7r  cos  4)q/Nq)  -|-  sin^(7rsin  (/)o/A^o) 
sin(7r/iVo) 


(3.27) 


The  parameter  q  in  (3.27)  is  At  independent.  In  other  words,  q  does  not  need  to  be  re¬ 
calculated  if  the  time  step  is  changed.  The  reason  we  still  use  the  least  squares  solution 
is  that,  for  the  3D  case,  the  zero  dispersion  angle  is  not  obvious  and  may  still  have  to  be 
numerically  calculated.  Substituting  (pQ  =  22.5°  into  (3.27),  the  difference  of  q  between  the 
analytical  solution  (3.27)  and  the  reference  solution  (3.26)  with  large  I  [\q{analytical)  — 
q{reference)W  is  plotted  in  Fig.  3.13.  As  shown,  for  practical  mesh  resolution,  the  reference 


solution  obtained  using  the  least  squares  method  is  a  good  approximation  to  the  analytical 
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solution. 

Now,  let  us  re-examine  the  least  squares  solution  of  (3.26).  Substituting  ai  and  b  into 
(3.26)  yields 


1 


q  = 


E 


Clj  — 


7  ZlLi  [sin^(vrcos  (pi/No)  -h  sin^(7rsin  (pi/No) 


(3.28) 


2=1 


sin(7r/No) 

which  reveals  the  similarity  between  the  least  squares  solution  and  the  analytical  solution  of 
(3.27).  As  I  approaches  infinity,  the  term  under  the  square  root  in  (3.28)  can  be  rewritten 


as 


1  ^  2  ^ 
lim  -  a(</>i)  =  lim  — 


TT  2 

=  -  a{<p)d(p 

TT  Jq 


(3.29) 


2=1  2=1 

which  is  the  average  of  a{<p)  for  (j)  G  [0,7r/2].  For  a  continuous  function  a((^),  using  the 
integral  mean  value  theorem,  we  obtain 

f-7r/2 


2  7^/2 
TT  Jo 


2 

a((/))#  =  -a{(t)o)  /  d(t)  =  a{(j)o) 

TT  Jo 


(3.30) 


for  a  constant  (j)o  £  [0,27r],  which  agrees  with  the  formulation  in  (3.27).  This  suggests  that 
the  zero  dispersion  angle  cpo  can  be  found  by  solving 


2  7^/^ 

a(0o)  =  -  /  a{4>)d(l)  (3.31) 

TT  Jo 

The  integral  in  (3.31)  still  needs  to  be  evaluated  numerically. 

Besides  the  local  dispersion  error  defined  in  (3.9),  the  2D  global  dispersion  error  is 
defined  as 


77'2D  _ 

^glb  - 


c2tt 


d(j) 


(3.32) 


In  Fig.  3.14,  we  plot  the  global  dispersion  error  when  the  design  mesh  resolution  of  NS22 


and  INS22  is  set  to  10.  As  shown,  near  the  design  mesh  resolution,  the  INS22  scheme 
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(/) 


Fig.  3.13.  Difference  of  q  between  the  analytical  and  reference  solutions. 


exhibits  an  improved  average  dispersion  distribution  compared  to  the  other  two  schemes. 
However,  the  low-dispersion  performance  of  the  NSFD  methods  is  still  relatively  narrow 
band. 


3.4.3.  Dispersion  Analysis  for  INS24 

The  optimal  q  values  for  the  INS24  scheme  can  be  derived  following  a  very  similar 
procedure  summarized  in  the  previous  section.  In  this  section,  the  3D  INS24  scheme  will 


be  considered  as  an  example.  The  derivation  starts  with  the  dispersion  relation  of  the  3D 
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Fig.  3.14.  Two-dimensional  global  dispersion  errors  of  the  (2,2)  schemes  (5  =  0.6,  iVo  = 
10, q=  1.00414816). 


INS24  scheme  which  can  be  written  as 


[||  sin(7r/iVo)  -  ^  sin(37r/lVo)]^  ■  2  f 


sin^^TT  S/Nq) 


sm 


N 


9  .  /  TT  COS  (f)  sin  9  k\  1  .  (  Sn  cos  cj)  sin  6  k 

-  sm  - — - sm 


1  2 


N  k  24: 


+ 


+ 


9  .  I  TT  sin  (j)  sin  9  k\  1  .  (  Stt  sin  (/>  sin  9  k 

N  k  “  N  k 


N  k 
2 


9  .  (  TT  COS  9  k\  1 


-  sm 


^  —  :rT  sm 


Nk  24 


37r  cos  9  k 
N  k 


(3.33) 


Similar  to  the  2D  case,  forcing  k  =  k  at  N  =  Nq  at  I  angles  along  i;i)(0°  ~  90°)  and  J  angles 
along  0(0°  ~  90°)  (note  that  the  dispersion  error  for  the  INS24  scheme  is  90°  periodic  with 
respect  to  cj)  and  9)  leads  to  an  overdetermined  linear  system 


Fp{q)=q^ -^=0,p=l,2,---IxJ 


(3.34) 
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Fig.  3.15.  Errors  of  q  for  the  3D  INS24  scheme. 


where 


(Xp  — 

+ 

+ 


9  /  vr  cos  sin  9j\  1  (  2>tt  cos  (pi  sin  0, 

-  sin  - — - - - 7  sin  ' 


No 


24 


No 


9  /  vr  sin  (pi  sin  9j\  1  f  2>tt  sin  (pi  sin  0, 

-  sin  - — - - - 7  sin  ' 


No 


24 


9  /  TT  COS  9j\  1  f  Stt  cos  9n 

-  sin  — — — - 7  sin  ' 


No  J  24 


No 


No 

1  2 


l2 


l2 


i  =  = 


h  = 


9  .  /  TT 


-  Sin 


1  .  /37r 


- 7  Sin 


No  24^“AiVo 


Then,  the  q  value  for  the  3D  INS24  can  be  solved  following  a  similar  procedure  summarized 
in  (3.24)-(3.26).  The  reference  solution  is  found  by  using  I  =  J  =  2,881.  The  errors  of  q 
versus  different  I  are  plotted  in  Fig.  3.15.  As  shown,  the  error  of  q  decreases  when  I  or  No 
increases.  Using  I  =  J  >  400,  the  error  of  q  is  smaller  than  10“®  for  A^o  =  10. 

The  correction  parameter  q  can  also  be  found  by  forcing  k  =  k  at  N  =  No  along  zero 
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dispersion  direction  {(j)o,6o)  of  the  INS24  scheme.  This  leads  to 


q  = 


+ 


9  / TT COS  4>o  sin  6q\  1  / Stt cos  4>o  sin  6q 

-  sin  (  - — -  —  —  sin  - — - 

Nq  y  24  I  No 


9  f  TT  sin  4>o  sin  0o  \  1  /  Svr  sin  (j)o  sin  ( 

- 

24 


+ 


8^“l 

V  No 

9  ■  1 
8®“  ' 

(  TT  COS  9o 

V  No  , 

9  .  / 

"  vr  \ 

No 


- 7  Sin 

24 


No 


1  2 


1/2 


Sin 


- 


1 

- sin 


\NoJ  24  \No 


Stt 


-1 


(3.35) 


However,  the  exact  values  of  (po  and  9o  are  not  obvious.  In  fact,  there  are  multiple  solutions 
for  the  zero  dispersion  direction.  For  convenience,  we  only  need  to  find  4>o  and  6o  along 
the  principal  axes,  for  example,  forcing  6o  =  0°  in  (3.35)  and  solving  for  (f>o,  or  vice  versa. 
Assuming  9o  =  0°  and  substituting  the  reference  q  value  with  I  =  J  =  2,881  into  (3.35), 
(jro  is  solved  and  plotted  in  Fig.  3.16.  As  shown,  (f)o  converges  to  27.07°  as  I  increases. 

Using  the  definition  in  (3.9),  the  local  dispersion  errors  of  Fang’s  algorithm  and  the 
NS24,  INS24  schemes  for  No  =  10  and  S  =  0.49  are  plotted  in  Figs.  3.17  and  3.18;  q  is 
calculated  using  I  =  J  =  401.  Only  the  dispersion  in  the  first  octant  is  plotted  due  to  the 
rotational  symmetry  of  the  curves.  As  shown,  both  NSFD  schemes  exhibit  lower  dispersion 
error  than  Fang’s  algorithm.  The  average  error  of  the  INS24  scheme  is  shifted  down  to 
zero  compared  to  the  NS24  scheme.  The  maximum  phase  differences  (peak-to-peak)  are 
the  same  for  all  three  methods. 


The  three-dimensional  global  dispersion  error  can  be  defined  as 

/*27r  rn 


^glb  - 


10  Jo 


sin  9d9d(j) 


(3.36) 


The  global  dispersion  errors  of  Fang’s,  NS24,  and  INS24  are  plotted  in  Fig.  3.19.  Notice 
that  the  INS24  scheme  exhibits  the  lowest  average  dispersion  error  near  the  design  mesh 
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I  (J=l) 

Fig.  3.16.  Zero-dispersion  angle  for  the  3D  INS24  scheme  (6*o  =  0°,5  =  0.49,  =  Nq  = 

10,  q{I  =  J  =  2,881)  =  1.00034998). 


resolution.  Moreover,  both  NSFD  methods  are  relative  narrow  band. 


3.4.4.  Numerical  Simulations 

In  this  section,  a  variety  of  numerical  simulations  will  be  performed  to  justify  the 
formulations  of  the  improved  nonstandard  finite-difference.  All  the  q  values  are  computed 
using  the  least  squares  solution.  All  the  open  problems  (radiation  and  scattering)  are 
truncated  with  an  8-layer  PML. 

In  the  first  simulation,  a  single-frequency  TM^  polarized  plane  wave  is  normally  inci¬ 
dent  upon  a  two-dimensional  dielectric  square  cylinder.  The  cylinder  has  a  side  width  of 
d  =  2Xq  (where  Aq  is  the  free-space  wavelength)  and  a  relative  permittivity  e,.  =  4.  An 
arithmetic  average  permittivity  (e^  =  2.5)  is  assigned  along  the  free-space-material  inter¬ 


face.  The  sk  and  q  values  are  different  for  different  materials.  For  example,  for  Nq  =  20, 
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Fig.  3.17.  Three-dimensional  local  dispersion  error  of  Fang’s  algorithm  (S  =  0.49,  = 

No  =  10). 


Fig.  3.18.  Three-dimensional  local  dispersion  error  of  the  NS24  and  INS24  schemes  (5  = 
0.49,  N  =  No  =  10,q  =  1.00034933). 


77 


Fig.  3.19.  Three-dimensional  global  dispersion  error  of  the  (2,4)  schemes  {S  =  0.49,  A^o  = 
10,  iV  =  6  ~  20,  g  =  1.00034933). 

it  can  be  solved  qi  =  1.00102838,  q2  =  1.00414816,  and  qs  =  1.00258172  for  the  mesh  reso¬ 
lution  in  the  free  space,  the  dielectric  material  and  the  interface,  respectively.  The  bistatic 
RCS  using  Nq  =  20  (in  free  space)  is  computed  by  Yee’s,  NS22,  and  INS22  schemes  and  is 
plotted  in  Fig.  3.21.  The  reference  solution  is  provided  by  Yee’s  algorithm  with  A^o  =  160 
(in  free  space). 

A  good  agreement  between  all  four  solutions  is  visually  indicated  except  for  the  angles 
close  to  180°.  To  compare  the  performance  of  the  three  algorithms,  the  following  L2  error 
is  defined  and  used. 


L2 


M 


(3.37) 


where  A®*™'  and  are  the  simulated  and  the  reference  amplitudes  of  the  linearized  RCS, 
respectively;  M  represents  the  number  of  the  observation  angles.  The  L2  errors  of  the  three 
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d 


Fig.  3.20.  Geometry  of  the  dielectric  square  cylinder  (d  =  2Ao). 


Fig.  3.21.  Bistatic  RCS  of  a  two-dimensional  dielectric  square  cylinder  {S  =  0.6,  Nq  =  20 
in  free  space,  M  =  361). 
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Fig.  3.22.  The  L2  error  of  the  linearized  RCS  of  the  dielectric  square  cylinder  versus  the 
free-space  mesh  resolution  (M  =  361). 

algorithms  using  different  design  mesh  resolutions  (A^o  =  10, 20,  and  40  in  free  space)  are 
plotted  in  Fig.  3.22.  As  shown,  the  INS22  scheme  has  the  lowest  L2  error.  However,  as 
the  design  mesh  resolution  increases,  the  performance  of  the  INS22  approaches  that  of  the 
NS22. 

For  the  second  simulation,  a  TE^  polarized  plane- wave  is  incident  and  scattered  by  a 
three-dimensional  PEC  square  plate  [2].  The  plate  has  a  side  width  d  =  5Ao  .  The  design 
mesh  resolution  is  selected  to  be  A^o  =  6.  The  reference  solution  is  predicted  by  Yee’s 
algorithm  with  No  =  40;  q  is  calculated  to  be  1.01556617.  The  yz-plane  bistatic  RCS  of 
the  plate  for  incident  angles  9i  =  0°  and  9i  =  30°  is  computed  by  Yee’s,  NS22,  and  INS22 
schemes  and  plotted  in  Figs.  3.24  and  3.25. 


A  good  agreement  between  all  four  solutions  is  indicated  visually.  The  L2  errors 
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Fig.  3.23.  Geometry  of  the  PEC  square  plate  {d  =  5Ao). 


Table  3.2.  Error  {S  =  0.55,  As  =  Ao/6,M  =  361) 


Scheme 

L2  Error  {6i  =  0°) 

L2  Error  {9i  =  30°) 

Yee 

24.6727 

14.6364 

NS22 

23.3063 

10.5742 

INS22 

24.1345 

6.0993 
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Observation  angle  0  (degrees) 

Fig.  3.24.  Bistatic  RCS  of  the  three-dimensional  PEC  square  plate  (yz- plane,  =  0°,  A^o  = 

6,5’  =  0.55,  M  =  361). 


-80  -60  -40  -20  0  20  40  60  80 


Observation  angle  0  (degrees) 

Fig.  3.25.  Bistatic  RCS  of  the  three-dimensional  PEC  square  plate  (yz-plane,  6i  =  30°,  Nq  = 
6,5’  =  0.55,  M  =  361). 
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Fig.  3.26.  Geometry  of  the  PEC  square  cavity. 

calculated  using  3.37  are  summarized  in  Table  3.2.  For  6i  =  0°  case,  the  NS22  scheme 
has  the  lowest  error  because  the  major  energy  in  the  scattering  pattern  is  along  0  =  0° 
where  NS22,  according  to  Fig.  3.12,  has  the  lowest  dispersion.  However,  for  the  9i  =  30° 
case,  the  INS22  scheme  has  the  lowest  error  since  the  major  energy  of  the  pattern  is  directed 
along  the  specular  angle  of  0  =  30°  where  the  INS22,  according  to  Fig.  3.12,  has  the  lowest 
dispersion  error. 

The  first  example  using  the  higher-order  (2,4)  schemes  simulates  the  mode  field  struc¬ 
tures  and  resonant  frequencies  of  a  two-dimensional  PEC  square  cavity  (TE^).  The  cavity 
has  a  side  width  d  =  1  m.  Fang’s,  NS24  and  INS24  schemes  with  As  =  1  m  are  used  for 
comparison.  For  the  extended  stencil,  image  theory  is  used  to  model  the  PEC  boundaries 
of  the  cavity  [61] .  The  cavity  is  excited  by  a  Gaussian  pulse  with  significant  energy  within 
0-0.5GHz.  After  updating  96,000  time  steps,  the  time  history  of  the  field  at  an  obser¬ 
vation  point  is  Fourier  transformed  to  the  frequency  domain.  The  absolute  errors  of  the 


resonant  frequencies  are  plotted  in  Figs.  3.27  and  3.28. 
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Fig.  3.27.  Absolute  error  of  the  resonant  frequencies  of  the  two-dimensional  PEC  square 
cavity  (As  =  0.1  m,  5  =  0.5,  /o  =  0.3  GHz). 


Fig.  3.28.  Absolute  error  of  the  resonant  frequencies  of  the  two-dimensional  PEC  square 
cavity  (As  =  0.1  m,  5  =  0.5,  /o  =  0.335  GHz). 
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In  Fig.  3.27,  the  design  frequency  is  set  to  be  0.3GHz  which  is  the  resonant  frequency 
of  the  TE20  mode;  q  for  this  case  is  solved  to  be  1.00026584.  As  shown,  the  NS24  scheme 
predicts  nearly  the  exact  solution.  The  error  of  the  INS24  scheme  is  higher  than  that  of 
the  NS24  scheme,  but  still  much  lower  than  that  of  Fang’s  algorithm.  The  explanation  is 
that,  for  the  TE20  mode,  the  wave  vector  is  directed  along  the  x  axis  where,  according  to 
Fig.  3.18,  the  NS24  scheme  has  no  dispersion  and  the  INS24  scheme  has  medium  dispersion 
error  compared  to  Fang’s  and  NS24  algorithms.  In  Fig.  3.28,  the  design  frequency  is  set 
to  be  0.335GHz  which  is  the  resonant  frequency  of  the  TE21  mode;  q  is  determined  to  be 
1.00041064.  In  this  case  the  INS24  scheme  has  the  lowest  error  because  the  wave  vector  is 
directed  along  26.6°  where,  according  to  Fig.  Fig.  3.18,  the  INS24  scheme  has  the  lowest 
dispersion  error. 

The  last  simulation  is  a  two-dimensional  four-element  array.  The  geometry  of  the  array 
is  shown  in  Fig.  3.29.  The  array  is  excited  by  four  single-frequency  soft  sources  {Ez)  which 
are  uniform  in  amplitude  and  in  phase.  Fang’s,  NS24  and  INS24  algorithms  with  Nq  =  6 
are  used  to  compute  the  array  factor.  Using  (3.26),  q  is  found  to  be  1.00196039.  In  Fig. 
3.30,  the  array  factor  power  patterns  versus  observation  angles  are  plotted  (note  the  90° 
pattern  is  periodic).  The  errors  calculated  using  (3.37)  are  summarized  in  Table  3.3.  It  is 
evident  that,  the  INS24  produces  the  most  accurate  result  based  on  the  observation  of  the 
pattern  or  the  comparison  of  the  L2  error. 

Table  3.3.  Error  {S  =  0.5,  As  =  Ao/6,  M  =  361) 


Scheme 

L?'  Error 

Fang 

0.07786753 

NS24 

0.03835138 

INS24 

0.02539144 
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Fig.  3.29.  Geometry  of  the  four-element  array  {d  =  5.5Ao). 


Observation  angle  (j)  (degrees) 


Fig.  3.30.  The  xy-plane  array  factor  power  pattern  of  the  four-element  array  (As  = 
Ao/6,5  =  0.5,  M  =  361). 
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In  this  section,  we  present  an  improved  NSFD  algorithm  by  introducing  an  extra 
degree  of  freedom  into  the  frequency-optimized  space  step  of  the  regular  NSFD.  The  INS22 
algorithm  exhibits  lower  dispersion  error  than  the  regular  NS22  and  Yee’s  algorithms  near 
the  design  mesh  resolution.  However,  this  low-dispersion  performance  is  narrow  band. 
The  correction  parameter  q  of  INS22  is  independent  of  the  time-step  which  avoids  the 
recalculation  when  the  time  step  changes.  The  computation  of  the  parameters  for  the  regular 
NSFD,  sk^Q  and  scuo,  is  simple.  To  apply  the  INS22  algorithm,  only  slight  modifications  are 
needed  to  Yee’s  algorithm.  The  increases  of  the  CPU  time  and  memory  consumption  are 
very  minimal.  The  INS24  scheme  shares  the  above  properties  of  the  INS22  scheme. 

Clearly,  the  INS  schemes  mitigate  the  average  dispersion  error  for  the  multi-dimensional 
simulations.  However,  the  INS  schemes  may  not  be  suitable  to  simulate  broad-band  signals 
since  the  low-dispersion  areas  are  limited  in  the  vicinity  of  the  the  design  frequency. 

3.5.  The  Generalized  NSFD  Schemes 

To  resolve  some  of  the  previously  stated  issues,  the  concept  of  the  INS  schemes  can  be 
generalized  by  introducing  more  degrees  of  freedom  into  the  space  stencils.  In  the  rest  of 
the  report,  these  schemes  are  referred  to  as  Generalized  NSFD  (GNS)  sehemes  [93].  The 
GNS24  scheme  will  be  used  as  an  example  to  demonstrate  the  analytical  procedure. 

3.5.1.  The  GNS24  Scheme 

It  has  been  reported  that  the  dispersion  error  of  Fang’s  FDTD(2,4)  algorithm  can  be 
drastically  reduced  by  properly  adjusting  the  coefficients  of  the  space  stencil  [58]  -  [65]. 
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However,  in  the  aforementioned  literature,  the  procedure  to  determine  the  optimal  coeffi¬ 
cients  becomes  cumbersome  for  complex  stencils.  A  more  straightforward  approach  is  to 
use  the  least  squares  method  introduced  in  the  previous  sections.  Now  the  GNS24  scheme 
will  be  presented  by  considering  a  2D  domain  with  the  square  mesh  (Ax  =  Ay  =  As).  The 
extension  to  three  dimension  is  straightforward. 

As  stated  in  Chapter  2,  the  fourth-order  of  accuracy  will  not  be  achieved  for  the  entire 
FDTD(2,4)  algorithm  due  to  the  second-order  error  in  the  time  domain.  Thus,  for  this  case, 
a  space  stencil  with  fourth-order  accuracy  may  not  be  necessary.  In  other  words,  a  space 
stencil  with  second-order  accuracy  is  sufficient.  This  allows  additional  degrees  of  freedom 
in  choosing  the  coefficients  of  the  space  stencil. 

In  the  2D  Cartesian  coordinate,  the  space  stencils  for  the  CNS24  scheme  can  be  ex¬ 
pressed  as  [60] 

rfj(A)  =  C‘T+1/2  -  /{-./2)  +C,(/;^3/2  -  /;-3/2)  _  ^ 

As 

where  Ci  and  C2  are  free  parameters  whose  values  will  be  determined  based  upon  the 
dispersion  analysis.  If  Ci  =  9/8  and  C2  =  —1/24,  (3.38)  reduces  to  the  regular  fourth-order 
stencil. 

From  Taylor  series  expansion,  it  can  be  shown  that  the  constraint  for  Ci  and  C2  to 
guarantee  second-order  accuracy  is  [60] ,  [93] 

Cl +2,02  =  I  or  Cl  =  1  -  3C2  (3.39) 

This  condition  allows  one  independent  free  parameter  which  can  be  utilized  to  reduce  the 
numerical  dispersion.  A  more  relaxed  constraint  would  be  [93] 


Cl  -|-  3C2  =  p  or  Cl  =  p  —  3C2 


(3.40) 


88 


where  p  is  an  extra  free  parameter.  This  condition  makes  the  stencil  in  (3.38)  nearly  second- 
order  accurate  (it  can  be  verified  that  p  is  very  close  to  unity)  but  provides  two  independent 
free  parameters  to  further  reduce  the  dispersion.  For  convenience,  the  schemes  governed  by 
(3.39)  and  (3.40)  are  referred  to  as  the  GNS24-1  and  GNS24-2  schemes,  respectively.  The 
values  of  the  free  parameters  Ci  and  C2  can  be  determined  following  a  similar  procedure 
demonstrated  in  Section  3.4.3. 

First,  the  dispersion  relation  of  the  GNS24  schemes  can  be  derived  and  written  as 
1  of  7tS  \  ^  /  TT  cos  (l)k\  ^  f  Stt  cos  (pk\ 

^  .  f7rsm4'k\  ^  .  f  Stt  sin  (p  k\'\ 

+  (3.41) 

where  S  is  the  Gourant  number,  N  is  the  number  of  cells  per  wavelength,  and  k/k  is  the 
ratio  of  the  numerical  wave  number  to  the  exact  wave  number.  Since  (3.39)  is  a  special  case 
of  (3.40),  we  only  present  (3.40)  as  an  example. 

Second,  forcing  k/k  =  1  aX  N  =  Nq  along  (()j  =  0°  ~  90°,  i  =  1  ~  /  in  (3.41)  leads  to 
an  overdetermined  nonlinear  system  which  can  be  expressed  as 

FiiCi,  C2)  =  (o-iGi  -|-  6iC*2)^  -|-  (qGi  -|-  (iiC*2)^  —  e  =  0,  i  =  1, 2,  •  •  •  ,  /  (3.42) 


(3.43) 
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where 


fi  =  aj  +  cf 
9i  =  bi  +  df 


hi  —  *2(^(libi  -\-  Cidi^ 


The  least  squares  solution  can  be  found  by  minimizing  the  square  error  defined  as 
follows  [92] 

I  I 

ll^(9)f  =  E  =  E  -  e)'  (3.44) 

i=\  i=l 

This  can  be  achieved  by  solving  the  following  equations 


G(Ci,C2)  = 


d\\F{q)\\^ 

_ 

dCi 

d\\F{q)\\^ 

dC2 

=  0 


(3.45) 


Gi(Ci,C2) 

G2(Ci,C2) 

Since  (3.45)  is  still  a  nonlinear  equation,  Newton’s  method  will  be  used  to  solve  for  Ci  and 
C2-  Assuming  C  =  [C'i,C'2]^  and  J  is  the  Jacobian  of  G,  the  procedure  is  summarized  as 
follows 


C”+i  =  C”  +  AC” 


(3.46) 


where 


AC”  =  -J-^G” 


(3.47) 


In  (3.47),  J  can  be  expressed  in  explicit  form  as 


dGi 

dGi 

a2||F(q)||2 

d^Fiq)p 

dCi 

dG2 

acA 

aCiCa 

dG2 

dG2 

d^F{q)\\^ 

d^F(q)p 

dGi 

dG2 

aGiC2 

90 


where 


d^F{q)r 

dCi^ 

dnF{q)r 

dCidC2 

dnF{q)r 

dCo^ 


i=l  L 

i=l 

i=l 


i=l 
I 


^y+F 

dCiJ  ^  * 

dFi  dFi 
dCi  dC2 


Fi 


dFi 

dC2 


+  Fi 


d^Fj 

dCi^ 

d^Fi 

dCiC2 

d^Fj 

dC2^ 


Notice  that  a  good  set  of  initial  values  for  the  iteration  is  Cq  =  [9/8, —1/24]',  i.e.,  the 
coefficients  of  the  regular  fourth-order  stencil.  The  coefficients  for  the  GNS24-1  scheme  can 
be  solved  similarly. 

Using  the  coefficients  with  large  /  (/  =  90,  001)  as  reference,  the  error  of  the  coefficients 
for  the  GNS24-1  and  GNS24-2  schemes  are  plotted  in  Figs.  3.31  and  3.32,  respectively. 
To  guarantee  the  errors  to  be  smaller  than  10“^,  I  =  401  will  be  used  to  calculate  the 
coefficients. 

As  an  example,  the  coefficients  for  Nq  =  10  and  Nq  =  20  are  computed  and  summarized 
in  Table  3.4. 


Table  3.4.  Gomputed  coefficients  for  the  GNS24  schemes  {S  =  0.5) 


No  =  10 

0^ 

II 

to 

0 

GNS24-1 

Cl  =  1.08867882 

C2  =  -0.02955961 

Cl  =  1.08630375 

C2  =  -0.02876792 

GNS24-2 

Cl  =  1.12864752 

C2  =  -0.04437762 

Cl  =  1.12592203 

C2  =  -0.04232427 

Using  these  coefficients,  the  local  dispersion  error  defined  in  (3.9)  for  each  of  the  GNS24-1 


and  GNS24-2  schemes  is  plotted  in  Fig.  3.33.  As  shown,  the  GNS24-1  scheme  exhibits  lower 
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Fig.  3.33.  Local  dispersion  error  of  the  GNS24  schemes  {N  =  No  =  10,  S  =  0.5). 

average  dispersion  error  than  the  FDTD(2,4)  scheme,  while  the  maximum  phase  error  (peak- 
to-peak)  is  greater.  However,  the  GNS24-2  scheme  presents  no  visible  dispersion  errors. 

The  global  dispersion  errors  defined  in  (3.32),  for  Nq  =  10  and  A^o  =  20,  are  plotted 
in  Figs.  3.34  and  3.35,  respectively.  It  is  evident  that  the  GNS24-1  scheme  generates 
smaller  global  errors  for  all  the  mesh  resolutions.  In  other  words,  the  GNS24-2  is  a  broad¬ 
band  scheme  in  terms  of  dispersion  error.  However,  compared  to  the  INS24  scheme,  the 
GNS24-2  scheme  is  a  narrow-band  scheme  whose  dispersion  error  is  negligible  at  =  Nq. 
The  dispersion-error  null  of  the  GNS24-2  scheme  is  precisely  controlled  by  the  design  mesh 


resolution  Nq. 
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3.5.2.  Absorbing  Boundary  Condition 

To  perform  the  simulations  in  an  open  region,  appropriate  absorbing  boundary  con¬ 
ditions  (ABCs)  must  be  applied  to  truncate  the  FDTD  domain.  One  of  the  most  popular 
ABCs  is  the  Perfectly  Matched  Layer  (PML)  proposed  by  Berenger  [8]  -  [9].  In  this  sec¬ 
tion,  the  effectiveness  of  Berenger’s  PML  for  the  GNS24  scheme  will  be  verified  following  a 
procedure  presented  in  [94] ,  [95] . 

A  two-dimensional  TE^  FDTD  test  domain  Q,t  is  truncated  by  Berenger’s  PML  with 
various  number  of  layers.  A  uniform  square  mesh  with  Ax  =  Ay  =  0.015m  is  used.  There 
are  100  cells  in  the  x  direction,  and  50  cells  in  the  y  direction.  The  Courant  number 
is  set  to  5  =  0.5.  A  benchmark  domain  with  the  same  constitutive  parameters,  is 
used  to  simulate  the  wave  propagating  in  an  unbounded  space.  The  benchmark  domain  is 
large  enough  so  that  there  is  no  significant  reflection  from  the  boundary  during  the  time  of 
interest.  The  test  and  benchmark  domains  are  shown  in  Fig.  3.36. 

To  excite  the  test  and  the  benchmark  domains,  a  hard-line  source  on  Hz  component 
is  placed  at  the  center  of  the  test  and  benchmark  domains.  The  source  should  have  a  very 
smooth  transition  to  zero  in  the  time  domain.  A  waveform  that  satisfies  this  characteristic 
was  suggested  by  [94]  and  repeated  here  as  follows 


Hz{at  center  of  domain)  =  < 


a[10  —  15cos(ti;i^)  -|-  6cos{uj2C)  —  cos(w3.^)]  ^  <  r 


^  >  r 


(3.48) 


where 


1  _q  ^  .  27rm 

“  =  'T  =  10  ,  ^  =  nAt,  ujrn  =  - ,  m  =  1,2,3 

o2i\J  7” 
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Fig.  3.36.  The  two-dimensional  test  and  benchmark  domains. 


The  time  history  and  the  frequency  spectrum  of  the  source  in  (3.48)  are  plotted  in  Figs. 
3.37  and  3.38,  respectively. 

The  errors  caused  by  the  tested  ABCs  are  obtained  by  subtracting  the  field  at  the 
observation  point  within  Qt  from  the  field  at  the  corresponding  point  in  Os.  The  local 
error  is  defined  as 


Eloc 


{i,j)  -H^ 

Hzmn.ax 


(3.49) 


where  H^max  is  the  maximum  magnitude  of  the  components  at  the  center  of  the  lower  air- 
PML  interface  and  HJ (i,  j)  are  magnetic  fields  in  the  test  and  benchmark  domains, 

respectively.  The  observation  points  are  selected  along  the  lower  air-PML  interface,  as 
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Fig.  3.37.  The  time  history  of  the  smooth  pulse. 


Frequency  (GHz) 


Fig.  3.38.  The  frequency  spectrum  of  the  smooth  pulse. 
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shown  in  Fig.  3.36.  The  global  error  is  defined  by 

Egib  =  (bi)]'  (3.50) 

*  j 

As  a  reference,  the  local  errors  of  Yee’s  algorithm  with  4,  8,  and  16  PML  layers  at 
t  =  lOOAt  are  plotted  in  Fig.  3.39.  For  Nq  =  20  (at  1  GHz),  the  coefficients  of  the  GNS24-2 
scheme  are  computed  to  be  Ci  =  1.12592203  and  C2  =  —0.04232427.  The  local  errors 
for  Fang’s  and  the  GNS24-2  schemes  are  plotted  in  Fig.  3.40.  As  shown,  the  local  errors 
generated  by  Fang’s  and  the  GNS24-2  scheme  are  about  the  same  level  as  that  of  Yee’s 
algorithm.  The  local  error  of  the  GNS24-2  scheme  is  slightly  greater  than  that  of  Fang’s 
algorithm. 

The  global  errors  versus  time  for  Yee’s  algorithm  are  plotted  in  Fig.  3.41.  The  global 
errors  for  Fang’s  and  the  GNS24-2  schemes  are  plotted  in  Fig.  3.42.  As  shown,  the  global 
errors  for  all  three  methods  with  4  or  8  PML  layers  are  about  the  same  level.  However,  for 
16  PML  layers,  the  global  error  of  Yee’s  algorithm  at  a  later  time  is  much  smaller  than  that 
generated  by  Fang’s  algorithm  and  the  GNS24-2  scheme.  This  probably  can  be  attributed 
to  the  larger  error  at  the  air-PML  interface,  which  forms  a  higher  noise  floor  for  the  two 
schemes  using  the  extended  stencils. 

From  the  above  numerical  experiments,  the  Berenger’s  PML  is  able  to  effectively  absorb 
the  outgoing  waves  simulated  by  Fang’s  and  the  GNS24  schemes.  For  all  the  simulations 
presented  in  this  report,  no  instability  was  observed. 

3.5.3.  Numerical  Simulations 

Two  numerical  simulations,  the  2D  PEG  square  and  the  four-element  array,  will  be 


repeated  for  the  GNS24  schemes.  The  geometry  and  the  parameters  used  in  the  simulation 


Fig.  3.39.  The  local  error  of  Yee’s  algorithm  (As  =  0.015m,  S  =  0.5). 


Fig.  3.40.  The  local  errors  of  Fang’s  algorithm  and  the  GNS24-2  scheme  (As  =  0.015m,  S 


I 


Fig.  3.41.  The  global  error  of  Yee’s  algorithm  (As  =  0.015m,  5  =  0.5). 
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Fig.  3.42.  The  global  errors  of  Fang’s  algorithm  and  the  GNS24-2  scheme  (As 
0.015m,  S'  =  0.5). 
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can  be  found  in  Section  3.4.4.  In  the  first  example,  the  absolute  errors  of  the  resonant 
frequencies  for  the  same  2D  cavity  are  plotted  in  Figs.  3.43  and  3.44.  The  design  frequencies 
of  Figs.  3.43  and  3.44  are  set  to  /o  =  0.3  GHz  and  /o  =  0.335  GHz,  respectively.  The 
computed  coefficients  for  the  GNS24-1  and  the  GNS24-2  schemes  are  summarized  in  Table 
3.5.  As  shown,  for  all  of  the  resonant  frequencies,  the  errors  generated  by  the  GNS24-1 
scheme  are  always  smaller  than  those  of  Fang’s  scheme.  The  GNS24-2  scheme  predicts 
precisely  the  resonant  frequency  for  each  case  when  the  design  frequency  is  set  to  the 
corresponding  resonant  frequency.  However,  the  errors  increase  rapidly  when  the  resonant 
frequency  deviates  from  the  design  frequency. 

Table  3.5.  Gomputed  coefficients  for  the  GNS24  schemes  (5  =  0.5,  As  =  0.1m) 


/o  =  0.3GHz 

/o  =  0.335GHz 

GNS24-1 

Cl  =  1.08871679 

C2  =  -0.02957226 

Cl  =  1.08951410 

C2  =  -0.02983803 

GNS24-2 

Cl  =  1.12864748 

C2  =  -0.04437761 

Cl  =  1.12953082 

C2  =  -0.04508125 

In  the  second  example,  the  array  factor  of  the  four-element  array  is  also  simulated 
using  the  GNS24-1  and  the  GNS24-2  schemes.  The  computed  power  patterns  are  plotted 
in  Fig.  3.45.  The  computed  coefficients  and  the  L2  errors  are  summarized  in  Table  3.6. 
Gompared  to  the  L2  errors  in  Table  3.6,  the  error  of  the  GNS24-1  is  larger  than  those  of  the 
NS24  and  INS24  schemes  but  smaller  than  that  of  Fang’s  scheme.  However,  the  GNS24-2 
scheme  predicts  the  most  accurate  results. 

The  numerical  simulations  in  this  section  clearly  justifies  the  conclusions  based  on  the 


dispersion  analysis.  The  extension  to  the  3D  GNS24  scheme  is  straightforward.  In  addition. 
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Fig.  3.43.  Absolute  error  of  the  resonant  frequencies  of  the  two-dimensional  PEC  square 
cavity  (As  =  0.1  m,  5  =  0.5,  /o  =  0.3  GHz). 


Fig.  3.44.  Absolute  error  of  the  resonant  frequencies  of  the  two-dimensional  PEC  square 
cavity  (As  =  0.1  m,  5  =  0.5,  /o  =  0.335  GHz). 
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Fig.  3.45.  The  xy-plane  array  factor  power  pattern  of  the  four-element  array  (As  = 
Ao/6,5  =  0.5,  M  =  361). 


Table  3.6.  Computed  coefficients  and  the  error  {S  =  0.5,  As  =  Ao/6,  M  =  361) 


Coefficients 

L2  Error 

GNS24-1 

Cl  =  1.09461960 

C2  =  -0.03153987 

0.03989388 

GNS24-2 

Cl  =  1.13483392 

C2  =  -0.04977567 

0.01759012 
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the  capability  of  the  GNS24  scheme  will  be  explored  in  the  following  two  section. 

3.5.4.  The  Angle-Optimized  GNS24  Scheme 

In  practice,  there  are  some  problems  in  which  the  dispersion  errors  in  certain  propa¬ 
gation  angles  are  dominant.  For  example,  the  simulation  of  an  highly  enlongated  cavity  or 
the  scattering  from  an  electrically  large  PEG  strip.  An  Angle-Optimized  FDTD  was  pro¬ 
posed  in  [57]  using  a  filtering  technique.  However,  to  achieve  the  same  goal  of  reducing  the 
dispersion  in  certain  angle  ranges,  a  more  straightforward  approach  is  to  use  the  GNS24-2 
scheme. 

Recall  that,  to  obtain  the  stencil  coefficients  of  the  GNS24  scheme,  k/k  is  forced  to 
unity  in  (3.41)  for  all  angles  (/>  =  1°  ~  90°.  If  we  force  k/k  =  1  only  along  the  dominant 
directions,  a  scheme  with  low  dispersion  within  the  desired  angle  range  can  be  expected. 
The  GNS24  scheme  presented  in  the  pervious  section  is  a  special  case  where  all  directions 
are  equally  important. 

Two  examples  are  presented  to  verify  the  above  argument.  In  Fig.  3.46,  the  local 
dispersion  errors  within  the  desired  angles  =  0°  ~  15°  are  plotted.  In  Fig.  3.47,  the 
dominant  angles  are  set  to  c/i  =  30°  ~  60°.  The  computed  coefficients  are  summarized 
in  Table  3.7.  As  shown,  the  regular  GNS24-2  scheme  {(/i  =  0°  ~  90°)  exhibits  an  more 

Table  3.7.  Gomputed  coefficients  of  the  angle-optimized  GNS24  scheme  (5  =  0.5,  A^o  =  10) 


Cl 

C2 

=  0°  ~  15° 

1.12862263 

-0.04436822 

(j)i  =  30°  ~  60° 

1.12867282 

-0.04438673 

=  0°  ~  90° 

1.12864733 

-0.04437755 
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evenly  distributed  dispersion.  For  the  Angle-Optimized  GNS24  (AOGNS24),  the  dispersion 
errors  along  the  desired  directions  are  much  lower  than  those  of  the  regular  GNS24-2  scheme. 
However,  for  the  angles  of  less  interests,  the  AOGNS24  suffers  from  much  greater  dispersion 
errors  than  the  regular  GNS24-2  scheme.  It  is  evident  that,  by  properly  adjusting  the 
coefficient  values  of  the  stencil,  the  low  dispersion  directions  of  the  AOGNS24  scheme  can 
be  flexibly  controlled.  Notice  that  the  variations  of  the  coefficient  values  for  different  angle 
ranges  are  subtle. 


3.5.5.  The  GNS24  Scheme  for  the  Rectangular  Mesh 

The  dispersion  analysis  conducted  so  far  are  all  using  a  square  mesh  (Ax  =  Ay  =  As). 
Now,  we  will  examine  the  performance  of  the  GNS24  scheme  for  the  rectangular  mesh 
(Ax  7^  Ay).  For  the  2D  case,  let  us  assume  the  ratio  of  the  space  increments  with  respect 
to  X  and  y  is  defined  as  R  =  Ax/ Ay  >  1.  Therefore,  the  mesh  resolutions  in  the  x  and  y 
directions  can  be  expressed  as  =  A/ Ax  and  Ny  =  A/ Ay  =  R  x  Nx,  respectively.  The 
Gourant  number  is  defined  based  upon  the  larger  space  step  (Ax)  as  S'  =  cAt/Ax. 

For  the  rectangular  mesh,  it  can  be  expected  that  the  coefficients  for  the  stencils  with 
respect  to  x  and  y  will  no  longer  be  identical.  Let  Cxi  and  Cx2  represent  the  coefficients 
for  the  stencil  with  respect  to  x  and  Cyi  and  Cy2  represent  the  coefficients  for  the  stencil 
with  respect  to  y.  The  2D  GNS24  dispersion  relation  for  the  rectangular  mesh  can  then  be 
written  as 


sin^ 


(^) 


1 


Ax^ 

1 

Ay2 


Cxi  sin 


Cyi  sin 


kxAx\ 
.  2  ) 
kyAy\ 


+  Cx2  sin 
-I-  Cy2  sin  I 


3kxAx\ 

—[ 

3kyAy\ 


(3.51) 


106 


where 


kx  =  kcos{4>),  ky  =  A;sin((^) 


Equation  (3.51)  can  be  rewritten  as 


ttS 


/  TT  COS  (j)  „  /  Svr  cos  (f)  k 

Cxi  sin  I  —  T  \  +  ^x2  sin 


Nx  k 


Nr  k 


^  .  (  TT  sin  (f)k\  ^  .  (  3tt  sin  cj)  k 


(3.52) 


The  procedure  to  solve  k/k  in  (3.52)  is  very  similar  to  the  square  mesh  case,  while  now  there 
are  totally  four  unknowns  Cxi-,  Cx2i  Cyi,  and  Cy2-  For  R  =  2,  Nxo  =  10,  and  S  =  0.25,  the 
coefficients  for  the  GNS24-2  scheme  are  computed  and  summarized  in  Table  3.8. 


Table  3.8.  Computed  GNS24-2  coefficients  for  the  rectangular  mesh  {R  =  2,Nxo  =  10,5  = 
0.25). 


Cxi 

1.05281421 

Cx2 

-0.01421612 

Cyl 

1.46453965 

Cy2 

-0.15900394 

The  local  dispersion  errors  of  the  Fang’s  algorithm  and  the  GNS24-2  scheme  are  plotted 
in  Figs.  3.48  and  3.49,  respectively.  It  is  evident  that  the  dispersion  curve  for  Fang’s 
algorithm  is  no  longer  90°  periodic.  Compared  to  the  curves  for  the  square  mesh  in  Fig. 
2.19,  a  high  dispersion  error  is  observed  in  the  vicinity  of  (j)  =  90°  due  to  the  different 
mesh  resolutions  in  the  x  and  y  directions.  However,  the  effects  of  the  rectangular  mesh 
is  automatically  compensated  by  the  GNS24-2  scheme.  A  very  low  and  more  isotropic 


dispersion  distribution  can  be  observed. 


Local  dispersion  error  (1-k7k)  -  2D 


Observation  angle  (degrees) 


Fig.  3.48.  Local  dispersion  error  of  Fang’s  scheme  (i?  =  2,  Nxo  =  10,  S  =  0.25). 


Fig.  3.49.  Local  dispersion  error  of  the  GNS24-2  schemes  (R  =  2,  A7o  =  10,5  =  0.25) 
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3.5.6.  More  Complex  Stencils 


The  procedure  summarized  in  the  previous  sections  can  be  applied,  but  not  limited  to, 
the  GNS24  scheme.  More  complex  stencils  can  be  derived  following  a  similar  procedure.  In 
this  section,  an  “Isotropic”  Finite  Difference  (IFD)  presented  in  [28]  will  be  considered  as 
an  example.  In  two  dimensions,  the  space  stencil  of  the  generalized  IFD  scheme  (GIFD) 
uses  six-point  stencils  which  can  be  expressed  as 


1,3) 


As 


+ 


^2(fi+l/2,j+l  +  fi+l/2,j-l  ~  fi-l/2,j+l  ~  fi-l/2,j-l) 


As 


(3.53) 


+ 


Dl{fi,j+l/2  -  fi,j-l/2) 

As 

D2{fi+l,j+l/2  +  /i-lj+1/2  “  fi+l,j-l/2  —  /j-l,j-l/2) 

As 


(3.54) 


where  Di  and  D2  are  undetermined  coefficients.  If  Di  =  11/12  and  D2  =  1/24,  (3.53)  and 
(3.54)  reduce  to  the  regular  IFD  scheme  in  [28].  The  2D  stencils  for  the  GIFD  stencils  are 
illustrated  in  Fig.  3.50. 

From  Taylor  series  expansion,  it  can  be  derived  that  the  constraint  for  (3.53)  and  (3.54) 
to  maintain  the  second-order  accuracy  is 


Di+2D2  =  1  OT  Di  =  l-  2D2 


(3.55) 


Similarly,  a  more  relaxed  constraint  is 

Di  +  2D2  =  q  or  Di  =  q  —  2D2  (3.56) 

where  q  is  an  extra  free  parameter.  We  refer  the  scheme  using  (3.55)  as  GIFD26-1  and  the 


scheme  using  (3.56)  as  GIFD26-2. 
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Fig.  3.50.  The  2D  general  six-point  stencils. 

The  dispersion  relation  of  the  general  (2,6)  scheme  can  be  derived  and  written  as 
1  .  2  f '^S\  r  f27rsm(pk\]  .  2  /vrcos (/> 

(wj  =  +  ™ 

^  ^  f  27rcos(^feM  9  f  TTsm6k\ 

+  +  (3.57) 

Following  a  similar  procedure,  I  equations  are  obtained  by  forcing  k/k  =  1  at  N  =  Nq 
along  01  =  0°  ~  •  ,  90°,  i  =  1,  2,  •  •  •  ,1.  This  leads  to  a  nonlinear  overdetermined  system 

which  can  be  expressed  as 

Fi{Di,  D2)  —  bi[Di  20*^2]^  -|-  di[Di  2ciD2]‘^  —  e  =  0,  i  =  1, 2,  •  •  •  ,  I  (3.58) 


no 


Equation  (3.58)  can  be  rewritten  as 

Fi{Di,D2)  =  fiDf  +  giD^  +  hiDiD2  -  e  =  0,i  =  1,2,  ■  ■  ■  ,I  (3.59) 

where 

=  bj  +  dj 

9i  =  4.{afbf  +  cjdj) 
hi  =  4{aibf  +  Cidf) 

The  resultant  overdetermined  system  in  (3.59)  can  be  solved  following  a  similar  procedure 
summarized  in  (3.44)-(3.48).  The  initial  values  are  chosen  to  be  the  coefficients  for  the 
regular  six-point  stencil,  i.e.,Do  =  [11/12,1/24].  The  computed  coefficients  using  I  =  401 
for  different  Nq  are  summarized  in  Table  3.9. 

Table  3.9.  Computed  coefficients  for  the  GIFD26  schemes  (5  =  0.5) 


No  =  10 

o 

CM 

II 

GIFD26-1 

Di  =  1.08464967 

D2  =  -0.04232484 

Di  =  1.08365977 

T>2  =  -0.04182988 

GIFD26-2 

Di  =  0.92662022 

D2  =  0.04292245 

Di  =  0.91913928 

T>2  =  0.04197646 

Using  the  computed  coefficients,  the  local  dispersion  errors  for  the  IFD26,  GIFD26-1, 
and  GIFD26-2  schemes  are  plotted  in  Fig.  3.51.  As  shown,  the  IFD26  scheme  does  exhibit 
an  “isotropic”  dispersion  distribution;  the  dispersion  variation  with  angles  is  negligible. 
However,  the  average  level  of  its  dispersion  is  still  very  high.  The  average  dispersion  error 
of  the  GIFD26-1  scheme  is  smaller  than  that  of  the  IFD26  scheme,  while  it  suffers  from 


Ill 


Fig.  3.51.  Local  dispersion  error  of  the  GIFD26  schemes  {N  =  Nq  =  10,  S  =  0.5). 

more  anisotropy.  As  a  comparison,  the  GIFD26-2  scheme  exhibits  no  visible  dispersion 
error  at  N  =  Nq. 

The  global  dispersion  errors  for  Nq  =  10  and  Nq  =  20  are  plotted  in  Figs.  3.52-3.53, 
respectively.  It  is  evident  that  the  GIFD26-1  algorithm  is  a  broad-band  scheme  in  terms  of 
the  average  dispersion  error.  However,  the  GIFD26-2  scheme  is  a  narrow-band  scheme.  At 
the  corresponding  mesh  resolution,  the  minimum  dispersion  errors  are  even  smaller  than 
those  of  the  GNS24-2  scheme  in  Figs.  3.34  and  3.35. 

3.5.7.  Stability  Gondition 

Another  important  issue  for  the  NSFD  methods  is  the  stability  condition  which  guar¬ 
antees  the  convergence  of  the  algorithm.  Gonventionally,  the  stability  condition  can  be 
obtained  by  performing  a  Von  Neuman  analysis  [88]  or  following  a  more  concise  approach 
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based  on  the  complex  frequency  analysis.  In  this  section,  the  procedure  of  [1]  will  be  used 
to  solve  the  stability  condition  for  the  NSFD  algorithms. 

Let  us  consider  the  2D  GNS24-2  scheme  as  an  example.  In  (3.51),  solving  for  lv  leads 


to 


c.=  ^sm-‘(0 


(3.60) 


where 


^  =  cAt 


1 


Cxi  sin 


kj:Ax 

2 

kyAy 


+  Cx2  sin 


SkxAx 


l2 


i2 


(3.61) 


To  guarantee  that  uj  in  (3.60)  has  a  real  value,  ^  in  (3.61)  must  satisfy  ^  <  1.  Since  Cxi  > 
0,Cx2  <  0  and  Cyi  >  0,Cy2  <  0,  the  maximum  value  of  ^  is  achieved  when  {kx^x)/2  = 
{kyAy) j2  =  71  j2.  Consequently, 


=  cAt\ 


\Cxl-Cx2)^  ^  {Cyl-Cy2)^ 


Rearranging  (3.62),  we  obtain 

At  < 


Ax^ 


1 


Ay2 


(3.62) 


^  /(C.1-C.2)"  ,  {Cyl-Cy2r 
^  Ax^  '  Ay‘^ 


(2D) 


(3.63) 


which  is  the  stability  condition  for  the  2D  GNS24-2  algorithm.  Following  a  similar  proce¬ 
dure,  the  ID  and  3D  stability  conditions  can  be  derived  and  expressed  as 

Ax 


At  < 


c(Ci  -  C2) 


(ID) 


(3.64) 


At  < 


(3D) 


,  KC^1-C^2?  I  {Cyl-Cy2?  , 

Ax^  ' 

For  the  cubic  mesh  case  (Ax  =  Ay  =  Az  =  As),  the  stability  conditions  reduce  to 


(3.65) 


S 


cAt  ^  1 

^  “  (Cl  -  C2) 


(ID) 


(3.66) 
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As  -  V2(Ci  -  C2) 


(2D) 


S=^<  — i - 

As  -  ^3(^1  -  C2) 


(3D) 


(3.67) 


(3.68) 


3.6.  Gauss’s  Laws 

Maxwell’s  equations  contain  Faraday’s  Law  in  (2.5),  Ampere’s  Law  in  (2.6)  and  Gauss’s 
Laws  for  electric  and  magnetic  fields  in  (2.7)  and  (2.8).  However,  most  FDTD  methods  only 
differentiate  the  former  two  equations  without  considering  the  two  Gauss’s  Laws.  Therefore, 
it  is  important  to  justify  that  any  FDTD  scheme  also  satisfies  Gauss’s  Laws.  In  source-free 
regions,  this  requires  that  no  artificial  electric  and  magnetic  charges  are  generated  inside 
the  domain.  In  Gartesian  coordinates,  it  has  been  verified  that  Yee’s  algorithm  satisfies 
Gauss’s  Law  for  rectangular  meshes  [1].  However,  it  remains  unclear  if  the  NSFD  methods 
meet  this  fundamental  requirement.  In  this  section.  Gauss’s  Law  for  the  electric  fields  will 
be  examined  for  the  NSFD  algorithms.  Gauss’s  Law  for  the  magnetic  fields  can  be  derived 
following  a  very  similar  procedure.  For  conciseness,  only  the  2D  (TE^)  case  is  considered. 
The  derivations  for  the  2D  (TM^)  and  the  3D  cases  are  similar. 

In  a  2D  free-space  region.  Gauss’s  Law  for  electric  fields  can  be  written  in  both  the 
differential  and  integral  forms  as 

Differential  form:  V  •  z3  =  0  (3.69) 

I  D  ■  dl  =  0 

c 


Integral  form: 


(3.70) 
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where  C  is  an  arbitrary  closed  contour.  The  time  derivatives  of  Gauss’s  Laws  in  (3.69)  and 
(3.70)  are  also  zero  since  all  the  fields  at  the  initial  time  are  zero.  Two  semi-discretized 
NSFD  update  equations  for  the  Ex  and  Ey  fields 

fO^-Ea;|j+l/2j  =  Dy^{Hz\iJ^i/2,j)  (3-71) 

^o^-E'i/ki+1/2  =  -Dx'^{Hz\ijj^i/2))  (3.72) 

will  be  used  in  the  derivation.  In  (3.71)  and  (3.72),  and  Dy^  represent  the  NSFD  space 
stencils  with  respect  to  x  and  y,  which  could  be  the  NS,  INS,  or  GNS  operators  introduced 
in  the  previous  sections. 

In  the  2D  Gartesian  coordinates,  the  rectangular  FDTD  grids  (Ax  7^  Ay)  are  shown 
in  Fig.  3.54.  For  clarity,  the  position  index  of  each  field  component  is  suppressed.  All  the 
involved  fields  are  numbered  in  sequence  according  to  the  superscripts.  Two  closed  contours, 
designated  as  Li  and  L2,  are  also  illustrated  in  Fig.  3.54.  Gauss’s  Law  in  differential  form 
is  defined  at  the  center  point  which  is  represented  by  the  empty  circle  in  Fig.  3.54  while 
Gauss’s  Law  in  integral  form  describes  the  total  charge  encompassed  by  the  Li  contour. 


3.6.1.  The  NS22  Scheme 


As  a  first  example,  let  us  consider  the  NS22  scheme  where  D”®  and  D”®  are  given  as 

(/i+l/2  “  fi-112) 


T~\ns22  E  _ 

ji  ~ 


-r\ns22  x  _ 


skx 

ifj+1/2  -  fj-1/2) 
sL, 


(3.73) 

(3.74) 


where  skx  and  sky  are  correction  parameters  which  possess  different  values  for  the  rect¬ 


angular  mesh.  Using  the  standard  second-order  central  difference,  the  time  derivative  of 
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Fig.  3.54.  Illustration  of  the  2D  rectangular  FDTD  grids  {TE^). 


where  E].,  E‘1  and  E^,  E^  are  used  to  evaluate  the  derivatives  in  the  x  and  y  directions, 
respectively.  Substituting  (3.71)  and  (3.72)  for  NS22  into  (3.75)  yields 


{Hi  -  Hi  -  Hi  +  Hi 


skyAx 


skxAy 


(3.76) 


Notice  that,  for  arbitrary  Tl^-fields,  Gauss’s  Law  will  not  be  satisfied  unless  the  follow¬ 
ing  constraint 

skx  sky 

Ax  Ay 


(3.77) 
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is  met.  It  is  evident  that  the  optimal  skx  and  sky  values  for  the  regular  NS22  scheme, 
defined  in  (3.4),  do  not  satisfy  this  constraint.  In  other  words,  the  regular  NS22  scheme 
does  not  satisfy  Gauss’s  Law  using  the  standard  differentiation.  If  we  were  to  enforce  the 
constraint  in  (3.77),  the  optimal  values  of  skx  and  sky^  defined  in  (3.4),  will  not  be  obtained. 
Consequently,  the  dispersion  performance  of  the  NS22  scheme  will  be  compromised. 

Naturally,  an  alternative  way  of  differentiating  (3.75)  is  to  use  the  NS22  operators  in 
(3.73)  and  (3.74),  which  are  self-consistent  to  the  scheme.  Following  a  similar  procedure, 
this  leads  to 


d 


skx 


sL, 


=  {Hl-Ht-Hl  +  R, 


1 


sky skx 


skxsk. 


=  0 


(3.78) 


That  is.  Gauss’s  Law  in  differential  form  is  satisfied. 

In  Cartesian  coordinate.  Gauss’s  Law  in  integral  form  can  be  evaluated  along  the  Li 
contour.  This  leads  to 


D-dl^  eo§-[{El  -  El)Ay  +  {EI  -  E^y)Ax]  (3.79) 

j  c 

Substituting  (3.71)  and  (3.72)  into  (3.79)  and  collecting  terms  yields 

4  £  Z5 .  -  ij;  +  ijj)  (^  -  ^)  (3.80) 

Clearly,  forcing  (3.80)  equal  to  zero  leads  to  the  same  constraint  as  in  (3.77).  Again,  if  Ax 
and  Ay  in  (3.79)  are  replaced  by  skx  and  sky^  respectively.  Gauss’s  Law  in  integral  form  is 


satisfied. 
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3.6.2.  The  GNS24  Scheme 


Next,  let  us  consider  the  2D  GNS24  scheme.  The  NS  and  INS  schemes  are  special 
cases  of  the  GNS24  scheme.  The  space  stencils  of  the  2D  GNS24  scheme  are  repeated  as 


D9ns2Aj.  ^ 

= 


Cxl{fi+l/2  —  fi-112)  +  Cx2{fi+3/2  —  fi-3/2) 


Ax 


Cyl{fj+l/2  -  fj-1/2)  +  Cy2{fj+3/2  -  fj-3/2) 
Ay 


(3.81) 

(3.82) 


Following  a  similar  procedure,  (3.69)  can  be  differentiated  using  the  standard  fourth-order 
stencil  and  written  as 


d 

d 


Ax 

(■EJ  -  E?)  -  A(EJ  -  E*„) 


27 
24\^y 


Ay 


(3.83) 


where  E^,  E'^,  E^,  E*  and  Ey^  Ey,  Ey,  Ey  are  used  to  evaluate  the  derivatives  in  the  x  and 
y  directions,  respectively.  Substituting  (3.71)  and  (3.72)  for  GNS24-2  into  (3.83)  yields 


[i27{Cy^-Cxi){Hl-H^-H^,+Hl) 

+  {27Cy2  +  Cxi){Hl-Hl^-Hl  +  Hl‘^) 

-  (27C,2  +  C,i)(77f -i7f  +  i7f) 

-  {Cy2  -  Cx2){h!  -  -  h!  +  Hl^)]  . 

Similar  to  the  NS22  case.  Gauss’s  Law  will  not  be  satisfied  unless  the  following  constraints 


Cxi  —  Cyi  —  Cl 
Cx2  =  Cy2  =  C2 


(3.84) 

(3.85) 


Cl  +  27^2  =  0 


(3.86) 
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are  met.  It  can  be  verified  that  the  optimal  GNS24  coefficients  in  Table  3.5.5  do  not  satisfy 
these  constraints.  Enforcing  these  extra  constraints  on  the  stencil  coefficients  compromises 
the  dispersion  performance. 

An  alternative  way  of  differentiating  (3.69)  is  to  use  the  GNS24  operators  defined  in 
(3.81)  and  (3.82).  This  leads  to  zero  divergence  and  can  be  expressed  as 


d 

\Cxi{El-El)  +  Cx2{El-Et) 

Ax 

d 

'Cyi{El-El)  +  Cy2{El-E^y) 

Ay 

=  0 


For  Gauss’s  Law  in  integral  form,  a  direct  integration  along  Li  leads  to 
D-dl  ^  {Cyi  -  -  Ht  -  Hi  +  Hi) 


m 


c 


+  Cy2{Hl-Hf-Hl  +  H, 

\  trQ  _i_  ijlO 


(3.87) 


(3.88) 


For  arbitrary  Hz  fields,  the  second  and  third  terms  in  (3.88)  will  never  vanish  unless  Cx2  = 
Cy2  =  0. 

However,  multiplying  both  sides  of  (3.88)  with  AxAy  and  reorganizing,  we  obtain 


{AxAy)^V-D 


d 


[Cxi{El  -  El) Ay  +  Cy^{El  -  eI)Ax] 


+  ^0 


'dt 

m 


^{El  -  El)3Ay  +  ^{El  -  EI)3Ax 


D  ■  dl  +  <b  D-dl 

L\  J  L2 


(3.89) 


Thus,  the  integral- form  of  Gauss’s  Law  can  be  approximated  using  the  weighted  sum  of  the 
closed  integrals  along  two  contours:  Li  and  L2.  Since  (3.89)  is  derived  from  Gauss’s  Law 
in  differential  form,  the  same  statements  can  be  made  for  Gauss’s  Law  in  integral  form. 


120 


That  is,  the  standard  formulation  does  not  satisfy  Gauss’s  Law  unless  extra  constraints 
are  enforced.  However,  using  the  nonstandard  formulation.  Gauss’s  Law  is  automatically 
satisfied. 

To  further  evaluate  Gauss’s  Laws,  a  simple  numerical  simulation  is  performed,  as  shown 
in  Fig.  3.55.  A  free-space  region  (TE^),  excited  with  a  soft  Gaussian  pulse  {ridecay  =  10, 
no  =  "iudecay  [1])  at  the  domain  center  {is,js),  is  simulated  using  both  the  NS22  and  the 
GNS24-2  schemes.  The  design  frequency  for  each  of  the  scheme  is  set  to  fo  =  300  MHz. 
The  simulation  domain  is  large  enough  so  that  there  is  no  reflection  from  the  edges  within 
the  simulation  duration  of  interest.  Rectangular  cells  with  Ax  =  0.1m  and  Ay  =  0.05m  are 
used.  The  Gourant  number  is  set  to  S'  =  0.25.  Electric  charge  densities  at  a  point  away 
from  the  source  {ig  +  3,  +  3)  are  computed  and  normalized  with  Dx{is,js)-  The  normal¬ 

ized  electric  charge  densities  are  plotted  in  Fig.  3.56  using  the  standard  differentiation  of 
(3.75)  and  (3.83)  in  Fig.  3.57  using  the  nonstandard  differentiation  of  (3.78)  and  (3.87). 
It  is  evident  that  artificial  charges  are  observed  in  the  source-free  region  if  the  standard 
formulations  are  used.  However,  the  nonstandard  formulations  reduce  the  artificial  charges 
by  a  factor  of  nearly  10^. 

Based  on  the  above  discussions,  the  following  conclusions  can  be  made  concerning 
Gauss’s  Laws  for  the  NSFD  algorithms.  In  the  source-free  region,  the  standard  formulations 
do  not  satisfy  Gauss’s  Law  unless  extra  constraints  on  the  stencil  coefficients  are  enforced. 
These  constraints  are  not  desirable  since  they  sacrifice  the  dispersion  performance  of  the 
NSFD  schemes.  Using  the  nonstandard  formulations.  Gauss’s  Law  (both  in  differential  and 
integral  forms)  is  automatically  satisfied. 
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3.7.  The  Material  Interface  Conditions 

One  of  the  greatest  barriers  for  the  schemes  using  the  extended  stencils  is  the  proper 
material  interface  treatments.  Even  if  there  is  no  staircase  error,  the  extended  stencils  fail 
to  accurately  predict  the  interaction  near  the  material  interfaces  [85].  Many  efforts  have 
been  made  to  minimize  the  errors  due  to  material  interfaces.  Among  them,  there  are  the 
hybrid  FDTD(2,2)/FDTD(2,4)  scheme  [72]-  [76],  the  one-sided  stencils  [77],  [78],  and  the 
derivative-matching  technique  [81].  Recently,  a  very  different  approach  was  presented  for 
Yee’s  algorithm  to  reduce  the  artifacts  due  to  the  material  discontinuities  [82]-  [84].  The 
approach  used  in  [82]-  [84]  is  similar  to  that  used  to  derive  the  NSFD  algorithms.  That  is, 
the  coefficients  in  the  space  stencils  are  allowed  to  be  adjusted.  By  forcing  the  numerical 
results  to  be  equal  to  the  exact  solutions  at  a  certain  design  frequency,  the  values  of  those 
free  coefficients  can  be  determined.  In  this  report,  the  basic  concept  of  [82]-  [84]  is  extended 
to  the  four-point  finite  difference  stencils. 
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Fig.  3.56.  Normalized  electric  charge  densities  using  the  standard  differentiation  of  (3.75) 
and  (3.83). 


Fig.  3.57.  Normalized  electric  charge  densities  using  the  nonstandard  differentiation  of 
(3.78)  and  (3.87). 
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Fig.  3.58.  One-dimensional  material  discontinuity. 


3.7.1.  Formulation 


In  a  one-dimensional  homogeneous  medium,  x-polarized  plane  waves  traveling  in 
directions  can  be  modeled  by  Maxwell’s  equations  in  the  form  of 


dE^  _  1  dHy 
dt  e  dz 
dHy  _  1  dE^ 
dt  li  dz 


(3.90) 

(3.91) 


Let  us  consider  that  a  dielectric  discontinuity,  formed  by  two  different  dielectric  media 
(ei  and  €2)  exists  at  2:  =  0,  as  shown  in  Fig.  3.58.  A  uniform  mesh,  with  the  space  step 
Az  is  used  everywhere.  Without  loss  of  generality,  an  E^  field  is  located  at  the  interface 
(2  =  0). 

The  exact  solutions,  for  the  normalized  E-  and  iL-fields  in  the  two  different  media,  are 


given  as  [2] 
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Medium  1: 

(3.92) 

(3.93) 

Medium  2: 

=  (1  -h  r)e“^'^22 

(3.94) 

^2  _  j^-.7fc22 

(3.95) 

V2 

where  r/i  and  t]2  are  the  wave  impedances  in  medium  1  and  2,  respectively;  ki  and  k2  are  the 
wave  numbers  in  different  media  at  a  design  frequency  /o;  and  T  is  the  reflection  coefficient 
from  medium  1  toward  medium  2  at  the  interface,  which  can  be  written  as 


r  = 


??2  -  m 


m  +  m 

In  a  homogeneous  medium,  the  regular  fourth-order  finite  difference  stencil  [15] 

/i-3/2  “  27/j_i/2  +  27/j_|_i/2  —  fi+3/2 


(3.96) 


m)  = 


Az 


(3.97) 


can  be  used  to  discretize  the  space  derivatives  in  (3.90)  and  (3.91).  However,  the  direct 
application  of  (3.97)  near  the  dielectric  interface  will  cause  significant  errors  [85].  In  [81],  it 
was  reported  that  these  errors  can  be  greatly  reduced  by  using  a  modified  stencil  near  the 
interface,  which  can  be  expressed  as 

C'l/i-3/2  +  C2fi-il2  +  C'3/i+l/2  +  C'4/i+3/2 


Dih)  = 


Az 


(3.98) 


where  Ci,  (72,  6*3,  and  (74  are  unknown  coefficients  which  need  to  be  optimized.  The  regular 
fourth-order  stencil  in  (3.97)  can  be  viewed  as  a  special  case  of  (3.98).  In  Fig.  3.58,  three 
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E'-fields  (at  2:  =  —Az,0,Az)  and  two  i^-fields  (at  z  =  —AzI^^AzjA)  will  be  influenced 
by  the  material  discontinuities  due  to  the  width  of  the  four-point  stencil.  Therefore,  the 
update  equations  for  the  E—  and  H—  fields  near  the  interface  can  be  modified  as 


(CI/+1  -  E\^) 

At 

= -1,0,1 

Ci  \  / 

(3.99) 

_  ^|n— 1/2^ 

=  --D{E\^),i  = -1/2, 1/2 
ta 

(3.100) 

At 

There  are  four  unknown  coefficients  for  each  of  the  five  update  equations  in  (3.99)  and 
(3.100).  In  [81],  the  method  to  find  the  unknown  coefficients  involves  the  use  of  fictitious 
points  and  solving  the  resultant  algebraic  equations.  Alternatively,  a  more  straightforward 
approach,  which  is  originally  proposed  for  Yee’s  algorithm  in  [82]  -  [84],  can  be  generalized 
to  optimize  the  general  four-point  stencil.  Unlike  the  treatment  for  Yee’s  algorithm,  two 
types  of  interface  conditions  can  be  derived  for  the  fourth-order  stencil. 

3. 7. 1.1.  Condition  1 

As  an  example,  let  us  consider  the  update  equation  for  the  £'-field  at  i  =  0.  Replacing 
the  E-  and  i7-fields  in  (3.99)  by  the  exact  solutions  in  (3.92)  and  (3.93)  at  a  design  frequency 
/o  and  equating  the  real  and  imaginary  parts,  we  obtain  two  equations  with  four  unknowns 
(Cl,  (72, 6*3,  (74).  By  forcing  Ci  =  1/24  and  C2  =  —27/24  (i.e.,  the  coefficients  for  the 
regular  fourth-order  stencil),  the  unknowns  reduce  to  two  (C3,  C4).  The  resultant  equations 
can  be  rewritten  in  the  matrix  form  as 


ail  «12 

C3 

61 

«21  «22 

C4 

62 

(3.101) 
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where 


24(1  +  r)  cos{k2Az/2)  24(1  +  F)  cos{3k2A.z/2) 

an  =  - ,ai2  =  - , 

m  m 

24(1  +  r)  sm{k2^z/2)  24(1  +  F)  sin(3A:2A2;/2) 

021  = - ,022  = - 


hi  = 


m  m 

(1  —  F)[27  cos(A:iA2;/2)  —  cos(3A:iA2;/2)] 


62  =  - 


m 

48(1  +  F)eAzsin(t(;At/2) 
At 


+ 


(1  —  F)[27  cos(A:iAz/2)  —  cos(3A:iA2;/2)] 


Vi 


(3.102) 


are  constants.  The  values  of  0^,0^  can  be  determined  by  solving  (3.101).  The  other  four 
sets  of  unknown  coefficients  for  i  =  ±l,±l/2  can  be  found  using  a  similar  approach.  In  the 
rest  of  the  report,  the  scheme  using  these  coefficients  is  referred  to  as  Condition  1. 


3.7. 1.2.  Condition  2 


Note  that  the  general  stencil  of  (3.98)  allows  a  total  of  four  free  coefficients.  However, 
in  Condition  1,  only  two  of  them  were  used.  To  take  advantage  of  the  extra  two  degrees 
of  freedom,  we  can  introduce  a  second  design  frequency  /g.  By  doing  this,  we  can  also 
expect  to  obtain  a  more  broad-band  scheme.  Again,  let  us  consider  the  update  equation  for 
the  if-field  at  i  =  0.  Substituting  the  exact  solution  into  (3.92)  and  (3.93)  for  two  design 
frequencies  /o,  /g  into  (3.99),  and  equating  the  real  and  imaginary  parts,  yields  the  matrix 
equation 


Oil 

012 

013 

Ol4 

021 

022 

023 

024 

a'li 

®12 

®13 

a'i4 

®21 

®22 

023 

®24 

Cl 

hi 

C2 

h2 

C3 

h'l 

C4 

.  ^2  . 

(3.103) 
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where 


an 

ai3 

021 

023 


bi 


(1  —  r)  cos(3A;iA2;/2)  (1  —  F)  cos(A:iA2;/2) 

- ,Oi2  =  - 

m  r/i 

(1  +  r)  cos{k2^z/2)  (1  +  r)  cos{3k2^z/2) 

- ,014  =  - 

m  m 

(1  —  r)  sin(3A:iAz/2)  (1  —  F)  sm{kiAz/2) 

- ,022  =  - 

m  Vi 


(1  +  F)  sm{k2^z/2) 


m 


024  — 


(1  +  F)  sin(3/c2A2:/2) 
V2 


0,62  =  - 


2Aze{l  +  F)  sin(a;At/2) 
At 


are  constants  corresponding  to  the  first  design  frequency  /q.  The  matrix  elements  j  and 
have  similar  expressions  as  atj  and  bij  except  that  Lv,ki,k2  are  replaced  by  to' ,k'i,k'2 
corresponding  to  the  second  design  frequency  /g.  The  values  of  Ci,  (72,  (7*3,  can  be  deter¬ 
mined  by  solving  (3.103).  The  other  four  sets  of  unknown  coefficients  can  be  determined 
similarly.  The  scheme  using  these  coefficients  is  referred  to  as  Condition  2. 


3.7.2.  Numerical  Simulations 


A  classical  problem  [77],  [81],  a  one-dimensional  partially-filled  dielectric  cavity  as 
shown  in  Fig.  3.59,  is  used  to  evaluate  the  aforementioned  interface  conditions.  In  Fig. 
3.59,  N  (even  number)  is  the  total  number  of  cells  used  to  mesh  the  cavity.  By  enforcing 
the  corresponding  boundary  conditions  at  the  PEC  and  dielectric  interfaces,  the  closed  form 
solutions  for  the  E-  and  Ff-fields  can  be  readily  derived  [2]  and  written  as 

sin(A:iz),  0  <  2;  <  1 

C  cos{k2z)  +Dsm{k2z),  1  <  z  <2 


E{z)  =  I 


(3.104) 


H{z)  =  { 


— cosikiz) 

AJIln  V  J-  / 


0<  z<l 


(3.105) 


— sin(A:22:)  -|-  sin(A:22;),  1  <  z  <  2 

V  Z,  /  I  JUJflQ  V  z  /5  - 
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.V 


PEC  Material  Interface  PEC 


Fig.  3.59.  One- dimensional  partially-filled  cavity. 

where 

^  n  ^  n  sin(/ci)cos(2/c2) 

(7  =  2  sinffei)  cos(k:2),  = - - - 

sin(A:2) 

Assuming  that  e^i  =  1  and  6^2  =  2.25,  one  of  the  resonant  frequencies  is  242.0138  MHz. 
This  frequency  will  be  used  in  all  of  the  simulations. 

The  regular  FDTD(2,4)  [15]  or  the  NS24  [91]  schemes  will  be  applied  to  the  update 
equations  in  the  homogeneous  regions.  Near  the  material  interface,  Condition  1  or  2  will 
be  enforced.  The  PEC  boundaries  are  modeled  using  image  theory.  The  only  differences 
between  the  regular  stencils  and  the  interface  conditions  are  the  coefficients.  Consequently, 
similar  fourth-order  stencils,  with  different  coefficients,  can  be  consistently  applied  through¬ 
out  the  domain.  The  coefficients  for  the  interface  conditions  will  be  calculated  once,  using 
(3.101)  or  (3.103),  before  the  time  loop  starts.  To  measure  the  accuracy,  the  error  is 
defined  as  [85] 

N 

11^112  =  .  AzY,{Ui-Ui)^ 

\  i=0 


(3.106) 


129 


Fig.  3.60.  The  snapshot  of  the  E-  and  i^-fields  for  Nq  =  80. 


where  Ui  represents  the  numerical  E-  and  //-fields  predicted  by  FDTD;  Ui  designates  the 
exact  solution  given  by  (3.104)  and  (3.105). 

In  the  first  example,  we  initially  consider  Condition  1.  Since  Condition  1  is  derived  for 
a  single  frequency,  NS24  will  be  applied  in  the  homogeneous  region  to  minimize  dispersion. 
In  the  simulation,  the  design  number  of  cells  and  the  Courant  number  are  set  to  be  Nq  =  80 
and  S  =  0.85,  respectively.  The  total  number  of  time  steps  for  different  N  is  2048(A^/24). 
Snapshots  of  the  E-  and  //-fields  for  Nq  =  80  are  plotted  in  Fig.  3.60.  It  is  evident  that 
the  NS24  scheme  with  Condition  1  predicts  almost  the  exact  solution.  The  error  of  the 
//-field  versus  N ,  with  and  without  Condition  1  near  the  material  interface,  are  plotted  in 
Fig.  3.61.  As  shown,  at  =  Nq  the  error  generated  by  NS24  scheme  with  Condition  1  is 
extremely  small.  However,  the  error  increases  rapidly  as  N  deviates  from  Nq. 
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As  a  second  example,  Condition  2  is  considered.  To  evaluate  the  broad-band  perfor¬ 
mance,  the  standard  FDTD(2,4)  is  used  throughout  the  homogeneous  region.  Indeed,  the 
update  coefficients  for  Condition  2  are  weak  functions  of  the  design  number  of  cells  Nq,Nq. 
In  other  words,  varying  A^o  and  Nq  will  not  alter  the  coefficients  significantly.  In  the  fol¬ 
lowing,  Nq  and  Nq  will  be  set  to  320  and  1  x  10®,  respectively.  To  demonstrate  the  order  of 
accuracy  in  the  entire  space  domain,  a  very  small  Courant  number  S  =  0.01  is  used.  The 
error  of  the  ili-field  versus  N,  with  and  without  Condition  2,  are  plotted  in  Fig.  3.62. 
As  shown,  due  to  the  error  generated  near  the  material  interface,  FDTD(2,4)  without  the 
interface  treatment  exhibits  only  second-order  accuracy.  However,  fourth-order  accuracy  is 
fully  restored  if  Condition  2  is  applied  near  the  interface. 

In  summary,  two  dielectric  interface  conditions  are  designed  for  the  four-point  space 
stencils.  Condition  1  is  almost  exact  at  the  design  frequency.  However,  this  ideal  perfor¬ 
mance  is  very  narrow  band.  Condition  2  generates  a  real  fourth-order  scheme  in  the  space 
domain.  However,  a  very  small  Courant  number  has  to  be  used  to  achieve  the  fourth-order. 
An  alternative  approach  is  to  use  a  fourth-order  discretization,  for  example  RK4  [85],  in 


the  time  domain. 


error  of  E  fields 


1o'  10^  10^ 

N 


Fig.  3.61.  error  of  the  Fi-fields  (Condition  1) 


N 


Fig.  3.62.  error  of  the  Fi-fields  (Condition  2) 


CHAPTER  4 


THE  EEEECTS  OE  PASSENGERS  ON  PED  MUTUAL  COUPLING 

An  important  EMI  issue  for  commercial  aircraft  is  the  possibility  of  interference  gen¬ 
erated  by  on-board  personal  electronic  devices  (PEDs),  such  as  cell  phones,  CD  players, 
laptop  computers,  etc.  “It  is  a  common  policy  of  all  commercial  airlines  to  prohibit  the  use 
of  PEDs  during  at  least  the  very  sensitive  phases  of  take-off  and  landing,  if  not  for  the  entire 
duration  of  the  flights.  The  policy  has  been  established  because  it  is  believed  that  radiated 
emissions  from  PEDs  can  interfere  with  the  aircraft  electronic  equipment,  i.e.,  jamming  in 
the  communication  systems”  [72],  [76]. 

In  Georgakopoulos,  et  al,  the  mutual  coupling  between  a  simulated  PED,  located  in 
the  cabin  area  of  a  simplified  scale  model  aircraft  fuselage,  and  an  antenna  mounted  on 
the  exterior  of  the  fuselage  was  predicted  using  the  EDTD  and  compared  with  measure¬ 
ments  [72],  [76].  This  work,  however,  was  performed  for  a  completely  empty  fuselage.  An 
airline  fuselage  is  of  course  highly  populated  with  various  structures:  metal  frames,  foam 
seat  cushions,  duct  work,  wires,  various  plastics,  insulating  and  decorative  materials,  etc. 
An  ideal  model  would  include  all  of  these  geometrical  details  and  material  properties,  but 
the  complexity  and  computational  resources  involved  in  such  an  undertaking  renders  it  im¬ 
practical.  Some  of  these  structures  can  be  neglected.  Surprisingly,  it  has  been  reported 
that  the  seat  cushions  and  even  their  metal  frames  have  a  negligible  impact  on  the  coupling 
of  such  antennas  [96].  Passengers,  on  the  other  hand,  are  a  class  of  structures  which  have 
a  high  probability  of  modifying  the  EMI  environment  found  within  the  perforated  metallic 
cavity  of  an  aircraft  fuselage  due  to  their  dispersive  and  lossy  constituent  materials.  Model¬ 


ing  the  enormous  complexity  of  the  human  body,  particularly  using  uniform  discretization 
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across  a  solution  space  on  the  scale  of  an  airliner,  is  also  prohibitive.  Instead,  the  low- 
order  approximate  passenger  modeled  here  is  represented  by  a  geometrically  simple  shape 
composed  of  a  single  bulk  material. 

In  this  chapter,  the  5-parameters  of  a  simulated  PED  located  within  a  scaled  fuselage¬ 
like  enclosure  and  an  antenna  mounted  on  the  exterior  of  this  enclosure  are  both  measured 
and  computed  using  the  FDTD  method  [1],  [6].  The  5-parameters  are  determined  for  the 
case  in  which  the  enclosure  is  populated  with  simplified  ’’passengers,”  and  compared  with 
those  for  the  case  of  the  empty  enclosure. 

4.1.  Geometry 

The  ’’Simplified  Fuselage”  is  an  aluminum  box  having  a  rectangular  cross  section  that 
would  enclose  the  fuselage  of  a  1:20  scaled  Boeing  757.  Although  individually  larger  than 
their  scaled  counterparts,  the  cabin  windows  have  approximately  the  same  aperture  area  per 
unit  length  as  those  of  a  scaled  757,  and  the  Simplified  Fuselage  has  a  large  aperture  at  one 
end  to  represent  the  cockpit  windows.  To  ease  fabrication  and  handling,  the  overall  length 
of  the  Simplified  Fuselage  is  75%  that  of  a  scaled  757.  Additional  dimensional  details  can 
be  found  in  [72],  [76].  The  image  in  Fig.  4.1  is  of  the  CAD  model  for  the  empty  Simplified 
Fuselage,  from  which  the  relative  locations  of  the  simulated  PED  and  the  external  antenna 
can  be  discerned. 

The  structure  of  the  human  body  is  extremely  complex,  and  consists  of  a  variety  of 
dispersive  and  lossy  materials  [1],  [11],  [97].  Modeling  all  the  details  of  the  human  body  is 
beyond  the  scope  of  this  work.  However,  approximately  70%  by  weight  of  the  human  body  is 
blood  which  has  salinity  close  to  0.9%  (grams/liter).  Therefore,  circular-cross-section  plastic 
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Fig.  4.1.  Geometry  of  the  empty  Simplified  Fuselage  (the  letters  denote  the  locations  of 
the  cabin  windows). 

tubes  containing  0.9%  salt  water  were  used  to  simulate  the  passengers.  As  illustrated  in 
Figs.  4.2  and  4.3,  there  are  6  passengers  on  every  row,  and  15  rows  (90  passengers  in  total) 
inside  the  fuselage.  The  distance  between  the  adjacent  rows  is  75  mm. 

4.2.  Measurements 

To  provide  experimental  data  with  which  to  compare  to  the  FDTD  simulations,  mea¬ 
surements  were  made  of  the  coupling  between  the  simulated  FED  and  the  external  antenna. 
These  measurements  were  made  with  the  Simplified  Fuselage  empty,  and  with  the  passen¬ 


gers  in  place. 
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Fig.  4.2.  Cross  section  of  the  Simplified  Fuselage  illustrating  the  geometrical  relationships 
between  the  simulated  FED,  the  exterior  antenna,  and  the  simulated  passengers. 


Fig.  4.3.  Geometry  of  the  Simplified  Fuselage  filled  with  90  passengers. 
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4.2.1.  Passengers 

The  low-order  approximation  to  scaled  passengers  consists  of  90  mm  long,  15.9  mm 
diameter  thin-walled  “propionate”  plastic  tubing  filled  with  0.9%  salt  water.  Top  and 
bottom  plug  caps  were  used  to  close  the  tubes.  A  small  hole  was  drilled  through  each  of 
the  top  caps  to  allow  the  air  displaced  as  the  cap  was  installed  to  escape. 

To  keep  the  passengers  in  their  designated  locations  and  upright,  an  expanded 
polystyrene  support  system  was  made.  In  a  sheet  of  1”  thick  polystyrene,  holes  were  cut 
according  to  the  FDTD  model  with  the  help  of  a  wood  template,  as  shown  in  Fig.  4.4. 
A  circular  brass  tube  with  its  edge  sharpened  was  used  to  cut  the  holes  into  which  the 
passengers  were  placed.  To  keep  the  passengers  from  falling  through,  a  second  layer  of 
expanded  polystyrene  was  placed  beneath  the  first.  Finally,  expanded  polystyrene  spacers 
were  used  beneath  the  support  system  to  place  the  passengers  at  the  designed  height  within 
the  Simplified  Fuselage.  The  complete  set  of  90  passengers  and  their  support  system  is 
shown  in  Fig.  4.5. 

The  permittivity  and  loss  tangent  of  the  saline  solution  was  measured  using  Hewlett- 
Packard’s  open-ended  coaxial  Dielectric  Probe  [98].  The  measurement  setup  is  shown  in 
Fig.  4.6.  The  measured  permittivity  and  effective  conductivity  of  the  0.9%  salt  water  at 
20°C'  are  plotted  in  Figs.  4.7  and  4.8,  respectively. 

4.2.2.  Measurement  Setup 

The  experimental  setup  is  essentially  identical  to  that  for  the  PED  measurements 
previously  reported  [72],  [76].  A  detailed  description  of  the  setup  and  measurements  follows. 

The  measurements  were  made  in  the  anechoic  chamber  to  preclude  the  possibility 
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Fig.  4.4.  A  top  view  photograph  of  the  expanded  polystyrene  passenger  support,  and  the 
wood  template  and  brass  tube  used  to  cut  the  holes. 


Fig.  4.5.  A  photograph  of  the  90  passengers  in  their  support  structures. 
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Fig.  4.6.  The  setup  for  the  measurements  of  the  permittivity  and  conductivity  of  the  0.9% 
salty  water. 


Frequency  (MHz) 

Fig.  4.7.  Measured  permittivity  of  the  0.9%  salt  water  at  20°C. 
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Fig.  4.8.  Measured  effective  conductivity  of  the  0.9%  salt  water  at  20°C'. 
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Fig.  4.9.  An  end-view  photograph  of  the  Simplified  Fuselage,  setup  in  the  anechoic  chamber 
with  the  instrumentation  beneath  it.  The  front  panel  had  not  yet  been  copper  taped  to  the 
fuselage,  revealing  the  simulated  passengers  and  their  expanded  polystyrene  supports. 

of  environmental  scattering.  The  Simplified  Fuselage  was  placed  near  the  center  of  the 
chamber,  and  was  supported  approximately  84  cm  (33  inch)  above  the  floor  with  blocks 
of  expanded  polystyrene.  The  only  significance  of  this  height  is  that  it  was  convenient  for 
placing  the  synthesizer  and  S-parameter  test  set  of  the  network  analyzer  beneath  the  model, 
thus  minimizing  the  lengths  of  the  RF  cables.  Carbon-loaded  foam  absorber  was  placed 
on  the  instrumentation  to  reduce  any  scattering  that  might  occur.  A  15  m  (50-foot)-long 
IEEE-488  cable  enabled  the  HP8510  processor  to  communicate  with  the  source  and  test  set, 
and  a  12  m  (40-foot)-long  IF  cable  plumbed  the  measured  results  to  the  HP8510  detector, 
located  outside  of  the  chamber.  A  photograph  of  the  setup  within  the  anechoic  chamber  is 


shown  in  Fig.  4.9. 
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One  of  the  RF  cables  passed  through  a  hole  in  the  fuselage  floor  and  into  the  hollow 
pedestal  on  which  the  monopole  (which  simulates  the  FED)  was  mounted.  A  51  cm  (20 
inches)  long  “Gold  Tipped”  Addams-Russell  precision  cable  was  used  here  to  capitalize  on 
its  interchangeable  interface.  A  female  APC-3.5  connector  was  used  to  measure  the  ’’Thru” 
standards  of  the  full  2-port  calibration.  The  female  connector  was  removed  and  replaced 
with  a  phase-matched  male  APC-3.5  connector  with  which  the  other  calibration  standards 
and  the  S-parameters  were  measured. 

A  1.5  m  (5-foot)-long  braided  cable  that  had  been  wrapped  with  ferrite-loaded  elas¬ 
tomeric  absorbing  material  was  used  to  connect  the  external  monopole  to  the  other  port  of 
the  test  set. 

4.2.3.  Calibration  and  Measurement  Parameters 

Extremely  rapid  amplitude  variations  occur  as  a  function  of  frequency  in  S'12  between 
the  PED  and  the  external  antenna  due  to  the  reverberations  within  the  fuselage.  However, 
the  maximum  number  of  frequency  points  that  can  be  measured  in  one  sweep  with  the 
network  analyzer  is  801.  This  is  insufficient  to  accurately  sample  these  frequency  response 
variations  over  the  desired  band  of  50  MHz  to  6  GHz.  To  increase  the  sampling  density,  six 
bands  of  801  frequency  points  were  measured.  The  first  band  ranged  from  50  MHz  to  1.0 
GHz,  the  remaining  bands  were  each  1.0  GHz  wide. 

In  order  to  measure  the  error-corrected  S-parameters  at  multiple  frequency  bands,  a  cal¬ 
ibration  had  to  be  performed  for  each  of  the  frequency  bands.  A  standard  open/short /load 
“full  2-port”  calibration  (with  the  Isolation  standard  omitted)  was  performed  at  each  fre¬ 


quency  band.  Since  the  HP8510  does  not  have  sufficient  memory  to  store  all  of  the  correc- 
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tion  coefficients,  the  six  801-frequency  calibration  sets  were  stored  on  3.5”  floppy  disks  and 
loaded  and  deleted  from  the  network  analyzer  memory  as  needed. 

The  other  instrumentation  parameters  of  note  are:  the  source  was  operated  in  “Step 
Frequency  Mode”  (the  oscillator  was  phase  locked  at  each  frequency) ,  the  source  power  was 
-|-10  dBm,  512  IF  averages  were  used  for  the  high-signal  calibration  standards  and  1,024 
IF  averages  were  used  for  the  low-signal  calibration  standards  and  for  the  measurements. 
Over  the  first  frequency  band  at  which  the  amplitude  of  S'12  is  lowest,  an  IF  averaging 
factor  of  2,048  was  used.  Good  measurement  practices  were  employed:  the  connectors  of 
the  DUT,  calibration  standards,  instruments,  and  RF  cables  were  cleaned;  RF  connections 
were  torqued  to  8  in-lb;  and  more  than  sufficient  time  was  allowed  to  elapse  with  the 
instrumentation  turned  on  and  with  the  DUT  in  place  to  ensure  thermal  stability 

4.2.4.  Measurement  Results 

The  measured  S'-parameters  of  the  fuselage  with  and  without  passengers  are  shown  in 
Figs.  4.10  -  4.12.  Because  the  system  is  reciprocal  (S'12  =  S21),  only  S12  is  plotted. 

From  Figs.  4.10  -  4.12,  it  is  evident  that  the  presence  of  the  passengers  has  a  negligible 
impact  on  Sn  (the  two  curves  are  indistinguishable),  which  is  the  reflection  coefficient  of 
the  exterior  antenna.  However,  the  passengers  did  have  a  significant  effect  on  S22  and 
S12.  As  shown,  S22  (the  interior  antenna)  and  S12  without  passengers  exhibit  much  more 
oscillatory  behavior  due  to  the  large  number  of  resonances  present  inside  the  fuselage. 
However,  S22  and  5i2  with  the  passengers  exhibit  much  smoother  distributions  since  many 
resonant  components  have  been  significantly  dampened  by  the  “human  bodies”  which  are 


highly  dielectric  and  lossy.  Due  to  this  phenomenon,  the  presence  of  the  passengers  has 
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Fig.  4.10.  Measured  S'!!  (the  exterior  antenna)  of  the  simplified  fuselage  with  and  without 
passengers. 


Fig.  4.11.  Measured  S22  (the  interior  antenna)  of  the  simplified  fuselage  with  and  without 
passengers. 
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Fig.  4.12.  Measured  512  of  the  simplified  fuselage  with  and  without  passengers. 

the  effect  of  slightly  reducing  the  possibility  of  aircraft  system  upset  by  FED  emissions  as 
compared  to  that  of  the  empty  fuselage. 

4.3.  Simulations 

The  mutual  coupling  between  the  internal  FED  antenna  and  the  antenna  mounted 
on  the  exterior  of  the  fuselage  was  simulated  using  the  EDTD  algorithm  [1],  [6].  The 
5-parameters  were  computed  following  the  procedure  described  in  [8]. 

Due  to  the  very  large  permittivity  values  of  the  salt  water,  a  finer  mesh  size  than  that 
used  for  the  case  of  the  empty  fuselage  was  necessary.  Therefore,  a  cell  size  of  Ax  =  Ay  = 
Az  =  2.5mm  (A/40  at  3  GHz  in  free  space  or  about  A/5  at  3  GHz  in  the  salt  water)  was 
used.  To  allow  the  excitation  pulse  to  decay  to  a  low  level,  the  simulation  was  executed  for 


60,000  time  steps. 
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From  the  numerical  experiments,  it  has  been  found  that  the  coupling  results  are 
not  very  sensitive  to  variations  in  the  passengers’  permittivity  and  conductivity.  To 
demonstrate  this,  5i2  values  were  computed  for  three  different  combinations,  e.g.:  the 
low-frequency  parameters  (e^  =  78.22,  cr  =  1.42S/m)  at  500  MHz,  the  high-frequency 
parameters  at  (e^  =  60.12,  ct  =  18.43S/m)  10,000  MHz,  and  the  average  parameters 
(er  =  70.35, (T  =  8.05S/m).The  5i2S  computed  using  these  three  different  sets  of  param¬ 
eters  are  plotted  for  the  frequency  band  of  1,500  MHz  to  2,500  MHz  in  Fig.  4.13  and  for 
the  frequency  band  of  4,500  MHz  to  5,500  MHz  in  Fig.  4.14. 

As  shown,  the  simulated  results  using  different  parameters  agree  to  each  other  very 
well  regardless  of  the  dispersive  characteristic  of  the  salt  water.  That  the  direct  path 
through  window  C  likely  accounts  for  the  majority  of  the  coupling  of  energy  from  the  FED 
to  the  exterior  antenna  probably  explains  this  insensitivity  in  Si  2  to  the  parameters  of 
the  passengers.  Based  upon  this  observation,  it  can  be  justified  that,  for  this  particular 
configuration,  a  complete  dispersive  modeling  of  the  salt  water  is  not  necessary.  Hence,  the 
average  permittivity  and  conductivity  values  (70.35  and  8.05  S/m)  were  used  for  the  rest 
of  this  paper. 

The  time  snapshots  of  the  fields  in  various  cross-section  planes  are  plotted  in  Fig. 
4.17  -  4.18.  The  brightness  of  the  plot  increases  with  the  magnitude  of  the  E^  field.  Clearly, 
the  E  fields  subject  significant  damping  inside  the  “passenger  bodies”.  Also  notice  the 
strong  E  fields  radiated  through  Window  C  in  Fig.  4.18. 

Moreover,  the  time  history  of  an  E^  field  in  the  free-space  region  inside  the  fuselage 
is  plotted  in  Fig.  4.19.  The  Ez  field  at  the  same  position  for  the  empty  case  is  also 


plotted  in  Fig.  4.19  as  a  comparison.  It  is  evident  that  magnitude  of  the  Ez  field  with 
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Fig.  4.13.  Simulated  S'12  using  different  electric  property  combinations  (1500  -  2500  MHz) 


Fig.  4.14.  Simulated  5'i2  using  different  electric  property  combinations  (4500  -  5500  MHz) 
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Passengers 


Fig.  4.15.  Logarithmic  plot  of  the  field  on  the  yz-plane  (front  view)  across  one  row  of 
passengers. 


Exterior  Antenna 


Fig.  4.16.  Logarithmic  plot  of  the  Ez  field  on  the  yz-plane  (front  view)  at  the  center  of 
Window  C. 
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Fig.  4.17.  Logarithmic  plot  of  the  field  on  the  xz-plane  (side  view)  across  one  column 
of  passengers. 


PED 


Fig.  4.18.  Logarithmic  plot  of  the  field  on  an  rcy- plane  (top  view)  at  the  fuselage  center 
in  the  y  direction. 
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Fig.  4.19.  Representative  time  history  of  the  field  inside  the  fuselage  with  and  without 
the  passengers. 

passenger  presented  decay  much  faster  than  that  without,  due  to  the  heavy  damping  inside 
the  passenger. 

The  measured  and  simulated  S'-parameters  with  the  passengers  are  plotted  in  Figs. 
4.20  -  4.22.  As  shown,  ^n,  which  is  essentially  just  the  reflection  coefficient  of  a  monopole, 
has  a  null  at  2500  MHz.  The  null  of  FDTD  simulated  result  is  shifted  to  a  lower  frequency 
by  about  100  MHz.  However,  S22,  which  is  the  reflection  coefficient  of  the  internal  antenna, 
exhibits  a  similar  trend  as  that  of  5ii  but  with  more  oscillations  due  to  the  fuselage’s  highly 
resonant  characteristic. 

Of  the  S'-parameters,  512  is  the  most  important  since  it  represents  the  energy  generated 
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by  the  FED  which  could  interfere  with  the  aircraft’s  communication  system.  To  better 
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Fig.  4.20.  Measured  and  simulated  Sn  of  the  simplified  fuselage  filled  with  passengers 
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Fig.  4.21.  Measured  and  simulated  S22  of  the  simplified  fuselage  filled  with  passengers 
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Fig.  4.22.  Measured  and  simulated  512  of  the  simplified  fuselage  filled  with  passengers. 

illustrate  the  variation  of  5i2,  Fig.  4.22  is  expanded  in  different  frequency  bands,  as  shown 
in  Figs.  4.23  -  4.26.  Very  good  agreement  is  observed  from  1,500  MHz  to  5,500  MHz.  The 
discrepancy  near  the  lower  frequencies  is  probably  due  to  the  extremely  low  coupling  level. 
However,  the  discrepancy  at  the  high  end  of  the  frequency  band  can  be  attributed  to  the 
large  dispersion  error  caused  by  very  discretization.  Moreover,  it  is  important  to  note  that 
the  maximum  coupling  level,  which  is  higher  than  -40  dB,  occurs  between  2.0  -  3.5  GHz. 
This  level  of  coupling  is  sufficiently  high  to  draw  the  attention  of  the  communication  system 
designers. 

4.4.  Conclusion 

In  summary,  the  effects  of  passengers  on  the  mutual  coupling  between  a  simulated  FED 


and  an  externally-mounted  antenna  were  investigated.  By  comparing  the  5-parameters  for 
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Fig.  4.23.  Measured  and  simulated  S12  of  the  Simplified  Fuselage  filled  with  passengers 
(1500  -  2500  MHz). 


Fig.  4.24.  Measured  and  simulated  S'12  of  the  Simplified  Fuselage  filled  with  passengers 
(2500  -  3500  MHz). 


Fig.  4.25.  Measured  and  simulated  S12  of  the  Simplified  Fuselage  filled  with  passengers 
(3500  -  4500  MHz). 


Fig.  4.26.  Measured  and  simulated  5i2  of  the  Simplified  Fuselage  filled  with  passengers 
(4500  -  5500  MHz). 
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the  empty  scale  model  fuselage  with  those  for  the  case  in  which  the  scale  model  is  popu¬ 
lated  with  90  bodies  which  simulated  human  passengers,  it  was  found  that  the  presence  of 
passengers  significantly  dampens  the  reverberations  that  occur  within  the  fuselage.  This 
dampening  effect  reduces  the  magnitude  of  the  coupling  over  most  of  the  frequencies  con¬ 
sidered,  and  by  up  to  approximately  15  dB.  The  reduced  threat  indicated  with  passengers 
present,  potentially  represents  a  relaxation  in  the  EMI  shielding  required  for  on-board  sys¬ 
tems  compared  to  that  indicated  by  predicted  and  measured  5-parameters  for  the  empty 
fuselage  case. 


CHAPTER  5 


SUMMARY,  CONCLUSIONS  AND  FUTURE  WORK 

The  major  objective  of  the  research  presented  in  this  project  was  to  develop  low- 
dispersion  finite-difference  time-domain  (FDTD)  methods.  In  particular,  such  methods  are 
of  great  interests  in  many  practical  applications  which  contain  electrically  large  objects 
and/or  require  very  long  simulation  time.  The  currently  standard  FDTD  method,  Yee’s 
algorithm,  is  only  second-order  accurate  in  both  the  space-  and  time-  domains  which  suffers 
from  serious  dispersion.  However,  Yee’s  algorithm  is  robust  and  very  easy  to  implement. 
Over  the  years,  numerous  practical  problems  have  been  successfully  simulated  using  Yee’s 
algorithm.  Therefore,  to  take  advantage  of  the  available  geometry  library,  the  ideal  method 
should  be  able  to  minimize  the  dispersion  errors  without  significant  modifications  based  on 
Yee’s  algorithm. 

5.1.  Summary  and  Conclusions 

In  Chapter  1,  the  issues  caused  by  the  electrically  large  problems  and  the  difficulties 
for  Yee’s  algorithm  were  described.  It  was  indicated  that  the  mesh  refinement  is  not  a 
very  efficient  way  of  reducing  dispersion.  Thus,  the  severe  dispersion  accumulation  of  Yee’s 
algorithm  greatly  restricts  the  electrical  size  of  the  structures  that  can  be  handled.  Based  on 
a  selective  survey  of  the  existing  literature,  the  existing  dispersion-reduction  FDTD  methods 
can  be  roughly  classified  into  two  categories:  the  standard  (SFD)  and  the  nonstandard 
(NSFD)  FDTD  methods.  The  SFD  methods  use  the  standard  central  difference  stencils 


which  are  derived  from  the  Taylor  series  expansion  by  minimizing  the  leading  truncation 
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error.  To  further  reduce  dispersion,  the  NSFD  methods  adjust  the  coefficients  of  the  SFD 
stencils  directly  based  on  the  dispersion  analysis. 

Chapter  2  introduced  the  fundamental  theories  of  the  SFD  methods.  The  governing 
equations  of  any  electromagnetic  phenomena,  Maxwell’s  equations  was  first  reviewed.  Then, 
the  procedure  to  derive  Yee’s  FDTD(2,2)  algorithm  was  demonstrated.  That  is,  the  E— 
and  H—  fields  are  staggered  in  both  the  time-  and  space  domains.  Then,  the  time-  and 
space-  derivatives  in  Maxwell’s  equations  are  directly  approximated  by  the  second-order 
central  difference  stencil.  To  improve  the  accuracy,  it  is  natural  to  consider  the  usage  of 
the  higher-order  finite  difference  stencils.  For  simplicity,  the  higher-order  stencils  were  only 
used  in  the  space  domain  which  constructs  the  FDTD(2,A^)  schemes.  The  approach  to 
derive  the  higher-order  stencil,  the  method  of  undetermined  coefficients,  was  demonstrated 
for  the  FDTD(2,4)  algorithm.  Then,  the  procedure  was  generalized  for  central  difference 
stencils  with  arbitrary  high  order.  It  was  indicated  that  the  highest  order-of-accuracy  that 
can  be  achieved  increases  with  the  number  of  points  involved  in  the  stencil.  Next,  the 
stability  conditions  of  the  SFD  methods  were  also  briefly  discussed.  Finally,  the  dispersion 
characteristics  of  the  SFD  methods  were  quantitatively  studied.  Based  on  the  dispersion 
analysis,  it  is  clear  that  the  dispersion  is  determined  by  the  space  and  time  increments, 
the  operating  frequency,  and  the  propagation  angles.  Using  the  higher-order  stencils  does 
help  to  reduce  the  dispersion  errors.  It  was  also  observed  that  the  local  dispersion  curve 
of  a  SFD  method  is  centered  about  its  own  average  value.  Moreover,  for  the  FDTD(2,N) 
scheme,  the  ideal  N-order  cannot  be  achieved  due  to  the  larger  errors  generated  in  the  time 
domain.  These  observations  can  be  exploited  to  further  reduce  the  dispersion. 


In  Chapter  3,  the  major  topic  of  this  project,  the  NSFD  methods,  were  presented. 
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The  discussion  started  with  the  fundamental  theory  of  the  NSFD  algorithm.  It  was  shown 
that  the  regular  NSFD  algorithm  was  constructed  by  simply  replacing  the  time  and  space 
increments  in  Yee’s  scheme  by  their  frequency-optimized  counterparts.  The  ID  NS22  scheme 
was  first  introduced  followed  by  a  detailed  dispersion  analysis.  A  few  numerical  simulations 
were  provided  to  verify  the  performance  of  the  ID  NS22  scheme.  It  was  observed  that 
the  NS22  scheme  is  a  narrow-band  scheme  in  terms  of  dispersion,  which  is  very  suitable 
for  single- frequency  simulations.  The  dispersion  null  can  be  flexibly  controlled  by  a  design 
mesh  resolution.  For  broad-band  signals,  the  ID  NS22  scheme  could  predict  results  no  worse 
than  Yee’s  algorithm  as  long  as  the  design  mesh  resolution  is  properly  selected.  Then,  the 
basic  concept  of  the  NSFD  was  extended  to  the  four-point  stencil  to  create  an  NS24  scheme. 
Dispersion  analysis  revealed  that  the  NS24  scheme  possesses  two  dispersion  nulls  of  which 
one  controlled  by  the  design  mesh  resolution  and  the  other  is  controlled  by  the  Courant 
number.  With  two  nulls,  an  ID  broad-band  scheme  can  be  designed  using  the  NS24  scheme. 

However,  these  dispersion  characteristics  are  only  limited  to  the  ID  case.  For  multi¬ 
dimensional  cases,  the  above  schemes  suffers  from  significant  anisotropy.  The  reason  is 
that  the  derivation  of  the  frequency-optimized  time  and  space  increments  did  not  consider 
the  wave  traveling  along  directions  other  than  the  principle  axes.  To  mitigate  the  average 
dispersion  error,  values  of  the  frequency-optimized  increments  might  me  slightly  adjusted. 
Following  this  idea,  an  improved  NSFD  (INS)  algorithm  was  presented  by  introducing  extra 
free  parameters  into  the  frequency-optimized  increments.  This  is  equivalent  to  shifting 
the  average  dispersion  values  of  the  regular  NSFD  towards  zero.  The  concept  of  the  INS 
scheme  was  applied  to  both  the  (2,2)  and  (2,4)  stencils  to  construct  the  INS22  and  INS24 


schemes.  The  values  of  these  free  parameters  can  be  written  in  closed  form.  However, 
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a  zero  dispersion  angle,  which  can  only  be  found  numerically,  must  be  specified.  A  more 
straightforward  approach  was  to  use  the  least  square  method  (LSM).  A  follow-up  dispersion 
analysis  indicated  that  the  INS  schemes  did  improve  the  dispersion  in  global  sense  near 
the  design  mesh  resolution.  However,  for  single-frequency,  the  local  dispersion  error  still 
depends  on  the  propagation  angles.  In  other  words,  for  certain  directions,  the  regular  NSFD 
methods  could  be  more  accurate  than  the  proposed  INS  schemes.  Numerical  simulations, 
including  cavity,  radiation,  and  scattering  problems,  were  presented  using  both  the  INS22 
and  INS24  schemes. 

The  INS  scheme  opened  a  door  for  further  dispersion  reduction  which  could  be  achieved 
by  introducing  more  free  parameters  into  the  stencil.  A  generalized  NSFD  (GNS)  was 
presented  in  this  spirit.  With  the  cost  of  the  order-of- accuracy  in  the  space  stencil,  the 
GNS  schemes  are  more  powerful  in  dispersion  minimization.  Two  GNS24  schemes,  one 
broad-band  and  one  narrow-band,  were  presented.  The  coefficients  of  GNS  schemes  were 
solved  using  the  LSM.  It  was  found  that  the  average  dispersion  errors  of  the  broad-band 
GNS  scheme  are  smaller  than  that  of  their  standard  counterparts.  However,  this  scheme 
suffers  from  more  anisotropy.  The  narrow-band  scheme  exhibits  extremely  low  dispersion 
at  the  design  mesh  resolution  regardless  of  the  propagation  angles.  However,  this  excellent 
performance  can  only  be  achieved  near  the  design  mesh  resolution.  In  order  to  conduct 
radiation/scattering  simulations,  Berenger’s  perfectly  matched  layer  (PML)  was  examined 
for  the  schemes  with  extended  four-point  stencils.  From  the  numerical  experiments,  it  can 
be  observed  that  Berenger’s  PML  is  able  to  effectively  absorb  the  outgoing  waves.  Two 
simulations,  a  cavity  and  a  four-element  array,  were  performed  using  the  GNS  schemes. 


which  verified  the  dispersion  analysis.  Moreover,  it  was  found  that  the  GNS  schemes  can  be 
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used  to  further  reduce  the  dispersion  within  certain  angle  ranges.  Also,  the  extra  anisotropy 
generated  in  a  rectangular-meshing  domain  can  be  automatically  compensated  by  the  GNS 
schemes.  Then,  the  successful  application  to  a  (2,6)  stencil  showed  that  the  proposed 
procedure  works  equally  well  for  more  complex  stencils. 

After  that,  the  stability  conditions  of  the  GNS24  schemes  were  derived  using  a  complex- 
frequency  analysis.  Gompared  to  their  standard  counterparts,  the  stability  limits  of  the 
NSFD  schemes  are  slightly  perturbed.  Furthermore,  Gauss’s  Laws  were  examined  for  the 
above  NSFD  methods.  It  was  proven  that  the  Gauss’s  Law  is  satisfied  if  a  cubic  mesh 
is  used.  Using  the  standard  finite  difference.  Gauss’s  Laws  are  conditionally  satisfied  for 
the  rectangular  mesh.  The  extra  constraints  will  compromise  the  dispersion  performance. 
However,  using  the  nonstandard  differentiation.  Gauss’s  Laws  are  automatically  satisfied. 
In  the  last  part  of  this  chapter,  the  error  generated  by  the  material  discontinuities  was 
investigated.  Two  conditions,  one  narrow-band  and  one  broad-band,  were  proposed  for  the 
ID  dielectric  interfaces.  For  the  narrow-band  condition,  the  simulation  showed  that  the 
error  generated  by  the  discontinuity  is  significantly  reduced  for  the  design  mesh  resolution. 
For  the  broad-band  condition,  the  fourth-order  accuracy  of  the  entire  scheme  is  fully  restored 
if  a  small  time-step  is  used. 

Ghapter  4  examined  an  important  EMG  problem  concerning  the  effects  of  passengers 
on  the  FED  mutual  coupling  in  a  simplified  Boeing  757  fuselage.  The  two  antennas  includes 
one  monopole  mounted  on  the  exterior  of  the  fuselage  and  the  other  in  the  interior.  The 
“bodies”  of  the  onboard  passengers,  which  are  highly  lossy  and  dispersive,  were  simplified 
and  modeled  using  salt  water  contained  in  sealed  plastic  tubes.  The  S-parameters  were 


simulated  using  Yee’s  algorithm  and  measured  in  the  Electro-Magnetic  Anechoic  Ghamber. 
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Both  the  simulated  and  the  measured  results  agree  to  each  other  very  well.  It  was  observed 
that  many  high-frequency  resonances  were  significantly  dampened  due  to  the  lossy  behav¬ 
ior  of  the  ’’human  bodies”.  The  presence  of  the  passengers  reduces  the  threats  from  the 
electromagnetic  interference  for  the  onboard  system. 

5.2.  Future  Work 

The  research  of  this  project  presented  a  successful  approach  to  minimize  the  dispersion 
without  significant  modification  of  Yee’s  algorithm.  Some  useful  recommendations  for  future 
developments  are  summarized  in  this  section. 

First  of  all,  the  NSFD  methods  presented  in  this  report  only  achieve  low  dispersion 
near  a  design  mesh  resolution.  Some  of  them,  such  as  the  GNS24-1  scheme,  exhibit  lower 
dispersion  over  the  entire  frequency  band.  However,  the  improvement  is  very  limited.  One 
of  the  most  straightforward  solutions  is  to  consider  even  more  complex  stencils,  such  as  the 
sixth-order  stencils.  Since  more  coefficients  (degrees  of  freedom)  are  involved,  the  approach 
presented  in  this  report  can  be  readily  extended  to  two,  three  or  even  more  design  mesh 
resolutions.  By  appropriately  choosing  these  design  mesh  resolutions,  a  controllable  broad¬ 
band  low-dispersion  scheme  can  be  expected. 

Second,  a  more  challenging  problem  is  to  derive  the  corresponding  material  interface 
conditions  for  the  complex  stencil.  The  dielectric  interface  conditions  proposed  in  this 
project  are  limited  to  only  ID  case.  It  is  still  necessary  to  derive  proper  interface  conditions 
for  multi-dimensional  cases,  where  the  reflection  and  transmission  mechanisms  are  more 
complicated.  Furthermore,  there  is  still  no  practical  PEC  interface  condition  for  extended 


stencils,  which  can  handle  complex  3D  problems  accurately  with  stability.  The  material 
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interface  treatments  for  extended  stencils  still  remain  one  of  the  most  challenging  problems 
in  the  FDTD  area. 

Finally,  the  NSFD  methods  discussed  so  far  only  considered  simple  materials.  For  com¬ 
plex  materials,  such  as  lossy,  dispersive,  and  nonlinear  materials,  the  approach  introduced 
in  this  project  cannot  be  directly  applied.  In  [42],  [43],  Cole  did  some  pioneering  work 
in  applying  NSFD  to  complex  material.  However,  they  are  still  not  practical.  Moreover, 
Berenger’s  PML  was  not  fully  justified  for  the  schemes  using  extended  stencils.  The  larger 
noise  basis  generated  at  the  air-PML  interfaces  could  cause  serious  problems  for  late  time. 
The  development  of  the  appropriate  ABCs  is  still  necessary. 
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